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Chapter 1 
INTRODUCTION 


Robert N. Colwell 


As we near the conclusion o£ our work under this Integrated Study, the 
question that becomes increasingly relevant is this: To what extent are the 

potential users , especially in the State of California , actively accepting 
and using this modem remote sensing technology? One can find highly favor- 
able but somewhat superficial answers to this question, including those 
documented with letters of praise written and signed by high officials in 
the very agencies that we seek most to serve.; One has difficulty, however, 
trying to find answers which are both highly favorable and highly rigorous. 

To paraphrase the assertion made in Chapter 6 of this report by our Social 
Sciences Group, the transfer of a technology (remote sensing), can be consider- 
ed completed only when that technology, being readily available in the market- 
place, becomes generally accepted practice by the user agency. And, when 
the chief officer of that agency, upon routinely assessing all available 
technology that might be brought to bear on a specific problem, decrees that 
the technology in question is the one that shall be used. To date there is vir- 
tually no instance of truly modern remote sensing technology which has been fully 
transferred to the managers oT" Ual if ornia ' s natural resourcbs. 

Amplifying on this statement, and the problems which technology transfer 
implies, our report states, in Chapter 6, that remote sensing technology will 
not transfer itself. Instead, there is commonly a 5-stage ’’adoptive process” 
by which this new technology is perceived, internalized and used: (1) awareness 

(2) interest (3) evaluation (4) trial, and (5) adoption. Since this process 
may take years, there frequently are problems in maintaining the necessary . 
momentum, especially when there usually are disruptions in personnel and 
support along the way that can undermine both credibility and morale. Even 
after remote sensing technology has been adopted, it still may take years before 
this new technology will begin to bear fruit. 

In addition, the following pitfalls are likely to be encountered: 

(1) Oversell results in the raising of user expectations to levels beyond 
what current capabilities can deliver; (2) Overkill results when the user 
is urged to use elaborate techniques of computer-assisted analysis even when 
the desired information could have been derived quite adequately through the 
use of simple, inexpensive, and more readily understood manual interpretation 
techniques; (3) Under training results when a novice who has just completed 
an ’’appreciation” course in remote sensing, is required to plunge directly 
into the demanding tasks that are involved in making operational use of modern 
remote sensing technology; (4) Under involvement results when the user agency, 
plagued by a lack of qualified and/or motivated personnel, turns over the bulk 
of the work to consultants or others who lack familiarity vjith the user 

agency's resource problems and information needs and perhaps even with the 
resource as well; (5) Spurious Evaluation results when the user agency. 



forced by those in higher authority or by others into a ’’rush to judgment,” 
produces premature, incomplete, incestually-validated and usually overly- 
optimistic appraisals; and C6) Misapplication results, in part, from the sheer 
glamour of the shiny new tool known as ’’remote sensing” and is perhaps best 
metaphorized by the saying, "give a small boy a hammer and he soon discovers 
that everything needs pounding.” 

In light of the foregoing, let us now return to the question, "to what 
extent are the potential users in the State of California actively accepting 
and using this modern remote sensing technology?” From the preceding dis- 
cussion, state and local resource managers require more than a few years to 
adapt to their specific needs a technology as complex as remote sensing. 
Therefore, our overall conclusion is that here in the State of California, 
solid progress is being made and at a rate as rapid as might have been anti- 
cipated or rightfully expected. Furthermore, we believe that we can attribute 
much of this progress to the very substantial effort that we have made, even 
from the inception of this "Integrated Study” nearly nine years ago, to involve 
the potential users (i.e. the resource managers) throughout the entire process 
of technology development and transfer. In relation to the ultimate technology 
transfer objective of our study, the importance of having involved the ultimate 
users can scarcely be overemphasized. Nor can this involvement be overemphasized 
as we seek to provide a truly meaningful answer to the question that comprises 
the central theme of this introductory chapter: to what extent are California's 
potential users actually using, or at least progressing toward, the use of 
this remote sensing technology? Highly specific answers to that multifaceted 
question will be found in Chapters 3 through 5 of this report. We trust that 
those generally favorable answers when viewed in the light of the chapter-by- 
chapter "backup information" in support of them will not be regarded as the 
overly optimistic products of spurious evaluations made in the ivory-towered 
or ivy-clad vacuums of our own academia. In this regard, let us emphasize 
that a two-fold major reason why we have had the services of a Social Sciences 
Group throughout this study is to help point the way toward the achievement 
of true technology acceptance and to help us differentiate true acceptance 
from fancied acceptance of the type that might result from equal parts of 
wishful thinking and unwarranted optimism by our team of remote sensing scien- 
tists. Typical of the warnings that we have been given in this matter, is the 
following excerpt from Chapter 6 of the present progress report: 

"Most potential users of remote sensing technology are simply unmoved 
by paper-and -pencil evaluation games. They recognize that externally-prepared 
benefit-cost ratios exclude many of the considerations most important to them. 
They see impacts on their own decision processes, job security, and organization- 
al behavior being overlooked and obscured behind voluminous but vacuous evalua- 
tive reports. The result for the technology developers is ofen an evaluative 
"boomerang effect” in which users perform their own sub^ctive assessments and 
conclude, for various reasons, that fruits from the technology are not worth 
their price." 

Each evaluation relative to technology acceptance that has appeared in our 
Progress Reports since the start of this "Integrated Study” has been made in 
the light of this admonition. 

Finally, we believe that the heavy concentration of our present efforts 


on the preparation of ''Procedural Manuals" is essential to the maintaining 
of momentum toward our longer term goal. That goal is to bring about the 
complete acceptance and intelligent day-to-day utilization of modern remote 
sensing technology by resource managers in California, as well as in other 
parts of the world. As exemplified by the material appearing in Chapters 
3, 4, and 5 of this, and the two preceding progress reports, procedural 
manuals relate specific ways in which remote sensing techniques can be 
used advantageously in the inventory and management of California's natural 
resources. Furthermore, in each instance the manuals seek to oescribe those 
means by highly specific, step-by-step procedures. In the development of each 
of these Procedural Manuals, care has been exercised by our investigators 
to ensure that they were responsive to the intended user's needs. This included 
care in writing the step-by-step procedures, not only in a manner that could be 
understood, but also in a manner that could not be misu nderstood. 

As we near the final iteration on most of these Procedural Manuals, 
we continue to solicit reactions to them both from California's resource 
managers and from all other readers of this Progress Report. 


Chapter 2 
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Chapter 2 

WATER SUPPLY SITOJES BY THE DAVIS CAMPUS GROUP 


Co-Investlgator ; Ralph Algazi 
Contributor; Mlnsoo Suk 


(Explanatory Note by R. N. Colwell, Principal Investigator of the Integrated 
Study that Is being performed under this grant) ; U 

Attention is invited to Special Study No. 3 In the Pec. *77 Progress 
Report entitled "Satellite land Use Acquisition and Applications to Hydrologic 
Planning." That report was prepared by the long-time Co-Investigator from 
the Davis Campus for this NASA-funded, integrated study. Dr. Ralph Algazi 
and by his associate Dr. Mlnsoo Suk. Furthermore, the material contained 
in that Special Study is highly compatible with the emphasis that has been 
given during the past several years by all of our study participants relative 
to the remote sensing of California’s water resources. Nevertheless, the 
material does not appear in that Progress Report as a complete chapter, even 
though the contributions from Algazi et al . have been so treated in the past. 
This change merely reflects the orderly transition that has been occurring 
during the past year (with NASA approval) as the Algazi team has been bringing 
to a close its studies as funded under this particular NASA grant, and trans- 
ferring its remote sensing-related activities to follow-on programs under 
separate funding. 

The change just described is regarded by all concerned parties, both within 
NASA and within the University of California, as a healthy transition. This 
is especially so because the newly directed efforts of the Algazi team are 
brought about by their having been designated as participants in an Applications 
Systems Verification Test relative to uses that can be made of remote sensing 
in the estimation of water supply . Tims ^thelr re-directed efforts represent 
a logical step toward bringing about the acceptance by user agencies of certain 
kinds of remote sensing technology that have been developed, at least in part, 
by the Algazi group during their several years of participation in this grant. 

It is for the reasons indicated above that (1) for the first time there 
is not a full chapter in this Progress Report detailing grant-funded activities 
of the Algazi group, and (2) there is nevertheless a substantial contribution 
from that group appearing in Chapter 7 of our semi-annual progress report of 
31 December, 1977. 
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INTRODUCTION 


The focus of the work that has been perfornied during the 
present reporting period, by the Remote Sensing Research Program 
personnel at the Berkeley Campus, has been on completion of the 
remote sensing-aided procedural manuals on water supply forecasts. 
This work is described in terms of stepwise procedures which might 
best be followed by managers of water resources who are interested 
in using remote sensing techniques as an aid for water resources 
inventory, development, and management. 

The procedural manuals which follow this section are en- 
titled: 


• 1 
i 

r 

1. 

Remote Sensing as a Aid in Determining the Areal Extent 
of Snow 

i 

' 1 

1 t 

2. 

Remote Sensing as an Aid in Determining the -Water 
Content of Snow 


3. 

Remote Sensing as an Aid in Watershed-Wide Estimation 
of Solar and Net Radiation 

n. 

4. 

Remote Sensing as an Aid in Watershed-Wide Estimation 
of Water Loss to the Atmosphere 

^■i 

t k 

5 . 

Outline for Remote Sensing as an Aid in the Inventory 
and Management of the Multiple Resource Complex for a 
Given Wildland Area 

TTie 

manual on snow areal extent estimation is based on 


s j manual analysis of Landsat imagery supported by aerial photography 

\\ and may be easily used by any resource manager. 

« - - Tlie manual on snow water content estimation is based on 

; I snow areal extent data and a physical model for conversion of 

snow-cover data to water content of snow. These conversion 
models may be refined by further research. 

Il f 

-i The manuals on solar and net radiation estimation and water 

loss to the atmonphere estimation are both based on a data bank 
; [ composed of Landsat digital data, NOAA-4 satellite thermal data, 

terrain variations, and climatic data obtained from ground meteoro- 
logical stations . 
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Since both o£ these manuals consider variations in climate 
and terrain conditions, they may be applied to watersheds with 
different geographical locations. However, it should be men- 
tioned that not all the climatic effects have been taken under 
consideration. 

These procedures for solar and net radiation estimation 
require sophisticated computer facilities and special software 
packages for Landsat and NOAA-4 satellite data processing. 

Therefore, their applications are limited to the availability 
of hardware facilities with required software programs. 

The mathematical models used for solar and net radiation 
are based on energy-balanced equations and the results are in 
agreement with the data collected on the test site. This 
remote sensing-aided procedure for daily estimation of radiation 
components is recommended for application under most climatic 
and geographic conditions. The level I mathematical model for 
water loss (evapotranspiration) estimation, which is fully 
described and applied to the study area, is based on a semi- 
ompirical model. More I'efined models (level II and III) need 
to be applied to the study area for more accurate performance. 

It should be mentioned that it is very difficult to make, and, to the 
best of our knowledge, presently there is no operational method 
available for evapotranspiration estimation from wildland areas. 

This is the first time that remotely-sensed spatial data have 
been applied to watershed-wide evapotranspiration estimation 
and the results have been adequately accurate to be used in 
water yield prediction models. 

Presently, with the cooperation of the U.S. Forest Service 
at Plumas National Forest, the Plumas County Planning Department, 
and the California Department of Water Resources, we are con- 
tinuing the use of remote sensing for inventory and management 
of the multiple resource complex within the Upper Middle Fork 
of the Feather River Watershed. 
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PROCEDURAL MANUAL FOR 

REMOTE SENSING AS AN AID IN WATERSHED-WIDE ESTIMATION 
OF SOLAR AND NET RADIATION 


INTRODUCTION 

Of the tremendous quantity of radiant energy emitted by the sun, 
only a very small portion is intercepted by the earth and its atmos- 
phere. This small portion is the ultimate source of almost all of 
the earth’s energy. This amount of incident radiation is known as 
the solar constant which is defined (U.S. Army Corp of Engineers, 

1956) as: "the intensity of solar radiation received on a unit area 
of a plane normal to the incident radiation at the outer limit of 
the earth’s atmosphere with the earth at its mean distance from the 
sun". The value of the solar constant is generally taken to be 
1.94 Langleys per minute which is based on the 1913 Smithsonian stan- 
dard scale. There has been evidence that this value is too low; a 
value of 2.00 ly/min. has been offered (Miller, 1955). 

The portion of the solar radiation which actually reaches the 
earth’s surface depends upon the transparency of the atmosphere. In 
the absence of clouds these amounts are quite constant. The atmos- 
pheric transmission coefficient varies from about 80 percent at the 
time of the winter solstice to about 85 percent at the time of the 
summer solstice (U.S. Army Gorp of Engineers,, 1956). Atmospheric 
transmission coefficients are based on the total insolation received 
at the earth's surface and include both the direct and diffused 
radiation. 

By far the largest variations (in the amount of solar radiation 
transmitted by the atmosphere) are caused by clouds. The transmitted 
radiation varies with type, height, density, and amount of clouds. 

The quantitative consideration will be discussed later. Since clouds 
are such a powerful controlling factor in radiative heat exchange, 
other minor factors such as hximidity of the air are often ignored. 

Similar to the clouds effects, the ground cover conditions such 
as a forest canopy exert a powerful controlling influence on net all- 
wave radiation exchange in hydrologic processes. However, the forest 
canopy has a different effect than that of clouds, particularly with 
respect to shortwave radiation. While both the clouds and trees 
restrict the transmission of iradiation, clouds are highly reflective, 
while the forest canopy absorbs much of the radiation. Consequently, 
the forest canopy tends to be warmed and in turn gives up a portion 
of its absorbed energy for evapotranspiration. 
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Objective 


The objective of this procedural manual is to describe a 
remote sensing-aided methodology developed by personnel of the 
Remote Sensing Research Program for accurate and location-specific 
estimation of radiation components. These components include net 
shortwave (solar) radiation, net longwave radiation, and net 
radiation on the earth's surface. The results of this study are 
designed to be readily incorporated into evapotranspiration and 
other hydrologic models for water yield forecasts (Khorram, et al . , 
1976) . Since the inputs to the radiation estimation models entail 
both climatic and physiographic conditions of the watershed, the 
effect of climatic variability has been implicitly accounted for. 

Definitions 

Total incoming radiation (D and I) is defined as the direct (D) 
and diffused or indirect (I) shortwave radiation reaching the 
ground through the atmosphere. Some of this incident radiation 
is reflected bach to> the atmosphere (R), which is called reflected 
radiation and the rest is absorbed by objects on the earth's sur- 
face. Part of the absorbed radiation is dissipated to the atmos- 
phere from the ground (G) as longwave radiation. In turn, a 
portion of this dissipated longwave radiation is absorbed by the 
clouds and atmospheric particles and part of it is returned to 
the ground (A). This returned longwave radiation from the atmos- 
phere to the ground (A) is called atmospheric radiation. The dif- 
ference between upgoing (G) and atmospheric (A) longwave radiation 
may be defined as the net longwave radiation. Net radiation (Q ), 
as it is shown in the following equation is the sum of the net ^ 
shortwave radiation (D+I-R) and net longwave radiation (-G+A) . 

= (D+I-R) + (A-G) (1) 

PROCEDURE 

Prior to radiation estimation, the study area must be defined 
and the watershed boundary delineated, All of the information must 
be in common register, in this case Landsat-based; therefore, the 
Landsat data should be geometrically corrected and the vratershed 
boundary transferred to a Landsat data format. The techniques for 
geometric correction and for watershed boundary determination are 
described in Appendix I. 

The procedure for estimating net radiation is composed of a 
series of mathematical models and techniques for providing the 
required input information to the models. The items of input 
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The resulting equation is: 

Q^j = (-0.24 + 0.158 X lO'^^T^) (0.2 - 0.8 in cal. Cni"^’inin"l (7) 


Step 3. Net Shortwave Radiation Computation 

The value o£ Q is estimated by Linacre (1967) 
by means of the modified Angstrom equation. 

% = ^ ¥> 

where Q. is the value of Q above the at- 

mosphere, c and d are empirical constants 
related to the atmospheric turbidity, and 
n and N are as previously defined, n is 
measured and values of N depend on the 
latitude and time of year, and are available 
in standard meteorological tables. The 
values of n may be gained through remote 
sensing techniques using NOAA satellite 
data; however, n is normally measured in 
the field. Table 1 shows the values for 
c and d. Linacre 's equation is not exact 
because it does not allow for the varying 
importance of cloudiness at different times 
of the day. 


The values of for a day or less (con- 
sidering slope and aspect) are estimated 
by the following equation (Frank and Lee, 
1966): 

I 1 o 

Q.= —;r (t_ - ti) sin 0' • sin 6 + — • cos 9' 

A Z ^ X IT 


cos S. 


(sin wt^ - sin wtp where 


I = solar constant 
o 

e = radius vector, ratio of the earth-sun distance 
at a particular time to its mean 

t^ = number of hours before solar noon (negative) 

= number of hours after solar noon (positive) 


( 8 ) 


(9) 
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information are physiographic data, climatic data, and som6 
standard meteorological constants. 

Step 1. Allwave Radiation Model 

The net radiation flux Qj^ = (1-a) 

where a is albedo, is the total flux of 

shortwave radiation from the sun and the 
sky, and is the net upward flux of 

longwave radiation . 

Step 2. Net Longwave Radiation Computation 

The net longwave radiation is computed as 
follows: 

" ^Ld ~ ^Lu 

-104 -2 -1 

Qj^d = 0.971 X 10 TjJ - 0.245 (cal. cm. min 

0.813 X lO'^^T^ (cal. Cm~^ * min“b 

where is the net flux of longwave 

radiation in cloudless conditions, 

is the downward flux of longwave radiation 
from the formula of Swinback, 1963, (Tj^) is 

the ambient temperature in degrees of Kelvin 
and is the upward flux of longwave 

radiation from a surface with unit emissivity 
(Linacre, 1968). 

Then 


V “ '?L E ^ I 


where there are n hours of actual sunshine in a 
daylength of N hours. The term b has been 
taken to be either 0. 10 (Penman, 1948), 0.20 
(Kramer, 1957), 0,24 (Impens, 1963), or 0.30 
(Fitzpatrick and Stern, 1965); a value of 0.20 
may be adopted as a compromise (Linacre, 1968) . 
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5 = solar declination 

w = angular velocity of the earth's rotation CIS 
degrees per hour) 

12 1 . 

— = — in radians 
ir w 

9 = arcsin (sin k • cos h • cos 9 + cos k • sin 9) 
in which 9 - terrestrial latitude 

k = slope inclination in degrees 

h = slope azimuth (aspect) in degrees from the north 


wt' = wt + a 


where 
a = arctan ( 


sin h -sin k 


cos h • cos 9 - cos h > sin k • sin 9 


•) 


( 10 ) 

( 11 ) 


Step 4. Albedo Calculation 

The reflection or albedo is the ratio of reflected to inci- 
dent radiation, usually expressed as a percentage. The values 
of albedo for the total range of solar radiation, and in some 
cases for the visible electromagnetic range only, are shown in 
Table 2. The figures given in that table illustrate what can 
be observed directly from an airplane. The sea appeal's darkest, 
relative to white strips of surf and sand dunes. Woods show 
up darker than fields, and snow-covered areas are light. One 
can see here how this affects the whole heat economy of the 
earth's surface and therefore all hydrologic processes. 

Albedo is influenced not only by the nature of a surface, but 
also by its moisture content at any time, i.e. wet surfaces 
appear darker than dry surfaces. 

The procedure for:making albedo calculations is based on a 
vegetation/ terrain analysis of the watershed. Once the ground 
cover is determined, an albedo index is assigned for each type 
of vegetation or ground cover (based on the published values 
of albedo appearing in Table 1) and the albedo map of the 
watershed is prepared. A stepwise procedure for vegetation/ 
terrain analysis is described in Appendix II, following this 
procedural manual . 
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Step 5. Net Radiation Determination 

Combining the equations for net shortwave and net long- 
wave radiation, the net radiation on the earth's surface can 
be estimated by the following equation: 

= (l~a) Qg (-0.245 +0.158 x 10"^°T^) (0.2 + 0.8 (12) 

The values of the various radiation components estimated 
by this method are site- specific and consider the effects of 
slope, aspect, latitude, cloud cover and elevation. 

Step 6. Preparation of Topographic Inputs to the Radiation 
% Models 

The topographic characteristics of the watershed are some 
of the major factors affecting the amount of solar radiation. 

These characteristics include elevation, slope and aspect. 

The elevation, slope, and aspect maps are constructed 
from the elevation contour data on U.S. Geological Survey 
topographic maps, 15’ series. The procedure involves the 
orientation of each Landsat picture element (pixel) in the 
watershed. Each pixel is defined in terms of its x, y, and 
z cooi’dinate dimensions. The x and y coordinates locate the 
center of the pixel while the z coordinate defines the 
average elevation. Each pixel can be "tilted” and/or rotated. 

The "tilt and rotation functions are considered to be the 
slope and aspect of the pixel, respectively. 

The detailed procedures for coordinate determination 
and topographic analysis are described in Appendix III. 

Step 7. Preparation of Climatic Input Variables to Radiation 
Models 

The collection and integration of climatic information 
are basic to the use of simulation models that are used in 
the estimation of solar and longwave radiation, and evapo- 
transpiration for a watershed. The manner in which the in- 
formation is collected, along with the type(s) of data used 
are dependent on the specific models used and the desired 
nature of the final output. Among all of the climatic 
factors, temperature plays a very essential role in radiation 
models. 

Lack of adequate ground meteorological stations, in 
general, creates a problem in satisfying climatic data input 
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requirements for radiation models. Therefore, we have chosen 
two alternative solutions: (1) computer-assisted interpolation 
of existing climatic data over the watershed; and (2) the use 
of the National Oceanic and Atmospheric Administration (NOAA-4) 
satellite Very High Resolution Radiometer (VHRR) data for 
the generation of temperature maps for the area. 

Collection and interpolation of ground meteorological 
stations data are composed of the following steps: 

a. Identification and selection of ground meteorological 
stations 

h. Registration of the data 

c. Recording and coding of the data 

d. Data interpolation for climatic parameter of interest 

e. Preparation of the data in a format to be used in 
simulation models 

f. Generation of final products 

The detailed description and the techniques for the 
above-mentioned steps are given in Appendix IV of this manual. 

Data processing techniques for the NOAA-4 satellite 
thermal (VHRR) data as well as a brief description of the 
VHRR scanner and data acquisition system and required sys- 
tem capabilities for NOAA-4 data processing are given in 
Appendix IV of this manual. NOAA-4 satellite data proces- 
sing entails the following steps: 

a. Data acquisition 

b. Location of the study area on the raw thermal data 

c. Geometric correction of the data 

d. Data enhancement and noise elimination 

e. Transformation of voltage values to temperature values 

f . Transformation of temperature data to the desired 
coordinate system 

g . Preparation of final temperature map of the area of 
interest 


Step 8. Preparation of Final Products 

After all the inputs (both constants and variables) to the 
radiation models have been prepared^ each of the radiation 
components may be estimated. These output products consist 
of shortwave radiation (net solar radiation) by using equation 
(8) , net longwave radiation by using equation (7) , and net 
radiation on the earth surfacfe by using equation (12). The 
output products are the areal distribution of parameters of 
interest for a given site and date and percentage distribution 
of data over the study ai'ea. 


Table 1. Published Values of the Factors c and d in Equation (9)- 



Source 

Location 

Latitude 

c 

d 

c+d 


CD 

(2) 

(3) 

C4) 

(5) 

(6) 



Black et al C19~4) 

Stockholm, Fairbanks 

59,65°N 

0.22 

0.52 

0.74 

Monteith (1966) 

Lerwick 

60°N 

0.23 

0.56 

0.79 

Penman (1948) 

Rot hams ted , U . K . 

52°N 

0.18 

0.55 

0.7-3 

Baier and Robertson (1963) 

Canada 

52°N 

0.25 

0.62 

0.87 

Black et al. (1954) 

Kew, U.K. 

51°N 

0.19 • 

0.57 

0.76 

Van Wijk (1963) 

Gembloux 

Sl'^N 

0.15 

0.54 

0.69 

It II II 

Versailles 

49°N 

0.23 

0.50 

0.73 

Tanner and Peltcn (1960) 

Wisconsin 

45°N 

0.18 

0.55 

0.73 

de Villeie (1965) 

El Aouina 

37°N 

0.28 

■ 0/45 

0.71 

de Vries (1958) 

Deniliquin 

36®S 

0.27 

0.54 

0.81 

Damagnez (1963) 

Tunisia 

35°N 

0.16 

0.59 

0.75 

Prescott (1940) 

Canberra 

35°S 

.0,25 

0.54 

0.79 

Black et al (1954) 

Dry Creek 

35°S 

0.30 

0.50 

0.80 

Page (1961) 

Capetwon 

34°S 

0.20 

0.59 

0.79 

Glover (1958) • 

Durban 

30^S 

0.25 

0.50 

0.76 

Yada (1959) 

Delhi 

29°N 

0.31 

0.46 

0.77 

Glover (1958) 

Pretoria 

26°S 

0.25 

0.50 

0.75 

»l fl 

Windhoek 

23°S 

0.26 

0.52 

0.78 

Page (1961) 

Tananarive 

19°S 

0.30 

0,48 

0.78 

Smith (1964) 

Jamaica 

18°N 

0.31 

0.49 

0.80 







Table 1, continued 


Source 

CD 

Location 

(2) 

Latitude 

(3) 

c 

C4) 

d 

(5).,^ 

c+d 

(6) 

Fitzpatrick (1967) 

Kimberly 

16°S 

0. 53 

0.43 

0.76 

Cackett (1964) 

Central .Africa 

15°S 

0.32 

0.47 

0.79 • 

Page (1961) 

Dakar 

15°N 

0.10 

0.70 

0.80 

Yada (1959) 

Madras 

13°N 

0.31 

0.49 

0,80 

Davies (1965) 

Kano 

12°N 

0.26 

0.54 

0.80 

Smith (1964) 

Trinidad 

11°N 

.0.27- 

0.49 . 

0.76 

Stanhill (1963) 

Benin City 

7°N 

0.26 

0,38 

0.64 

Davis (1965) 

Accra 

6®N 

0.30 

0.37 

0.67 

Black et a.l (1954) 

Batavia etc. 

6°S 

0.29 

0.29 

0.58 

Page (1961) 

Leopoldville 

4'^S 

0.21 

0.52 

0.73 

n II 

Singapore 

1°N 

0.21 

0.48 

0.69 

Glover (1958) 

Kabete, Kenya 


0.24 

0.59 

0.33 

Page (1961) 

Stanleyville 


0.28 

0.40 

0.68 

Rijks (1964) 

Kampala 

0° 

0.24 

0.46 

0.70 

Mean Values 

, 


0.245 

0.507 

0.752 

Hounam (1963) 

Australia 

12-43"S 

0..26 

0.50 


Davies (1965) 

West Africa 

5-15 N 

0.19 

0.60 


Black et al (1954) 

general 

- 

0.23 

0,48 


Page (1961) 

tropics 

- 

0.23 

0.52 


Turc (1961) 

general 

- 

0.18 

0.62 



Davies (1965) gave 0.28 and 0.33 for c and d respectively 
SOURCE: Linacre, 1968 


TABLE 2 


Stand 

Albedo Percent 

Source 

Fresh snow cover 

75-95 

Geiger 

Dense cloud cover 

60-95 

II 

Old snow cover 

40-70 

II 

Clean £irn snow 

50-65 

It 

Ice, sparse snow cover 

69 

Smithsonian 

table 

Clean glacier ice 

20-50 

II 

Light sand dunes, surf 

30-46 

Geiger 

Sandy soil 

15-40 

n 

Meadow and fields 

12-30 

n 

Meadow, low grass 

15-25^ 

Smithsonian 

table 

Field, plowed, dry 

20-25^ 

II 

Densely huilt-up area 

15-25 

Geiger 

Woods 

5-20 

M 

Dark cultivated soil 

7-10 

II 

Douglas -fir 

13-14 

Avery 

Pine 

14 

Brooks 

Conifers 

10-15 

Budyko 

Deciduous forest , fall 

15^ 

Krinov 

Deciduous forest, summer 

10^ 

II 

Coniferous forest, summer 

ga 

II 

Coniferous forest, winter 

3^ 

II 

Meadow dry grass 

10^ 

II 

Field Crops, ripe 

15^ 

II 

Spruce 

8-9 

Baumgartner 

Earth roads 

3^ 

Krinov 

Black top roads 

8^ 

Sewing Handbook 

»l It ‘ M 

9® 

Krinov 


Table 2,_ continued 


Stand 

Concrete road 
Buildings 

Mountain tops, bare 

Clay soil dry 

Clay soil wet 

Clouds, stratus overcast, 
0-500 feet thick 

Clouds, stratus overcast, 
500-1,000 feet thick 

Clouds, stratus overcast, 
1,000-2,000 feet thick 

I Clouds, dense, opaque 

I Clouds, dense, nearly 

I opaque 

i Clouds, thin 

Clouds , stratus , 600- 
1,600 feet thick 

i Clouds, stratocumulus 

overcast 

Clouds , altostratus , 
i occasional breaks 

1 Clouds , altostratus , 

! overcast 

Clouds, cirrostratus and 
altostratus overcast 

Clouds, cirrostratus 
overcast 

iWhole earth: 


Albedo Percent 

Source 

35® 

Krinov 

9® 

M 

24® 

If 

15® 

Sewing Handbook 

7.5® 

II 

5-63 

II 

31-75 

II 

59-84 

U 

55-78 

M 

44 

11 

36-40 

n 

78 

II 

56-81 

II 

17-36 

II 

39-59 

II 

49-64 

II 

44-50 

II 


From sjij;%S'urenients on the bright and dark portions of the moon, 
Danjon^*'' has calculated the albedo of the earth in the visible 
portion of the spectrum at 39 percent. Fritz^- has extended 
the calculation to include infrared and ultraviolet radiation, 
obtaining a total albedo of the earth of 35 percent. Baur 
and Philipps^ have calculated the total albedo to be 41.5 percent 


Visual albedo 

^Danjon, A., Ann. I'Obs. Strasbourg 3, No. 3, 1936, P. 139 

^Fritz, S., Bull. Amer. Meteorol . Soc., Vol. 29, 1948, p. 303; 

\^ol, 31, 1950, p. 251; Journ . Meteorol . , Vol. 5, 1949, p.277 

'^Baur, F., and Philipps, H., Gerl. Beitr. Geophys., Vol. 42, p,160 
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GEOMETRIC CORRECTION AND WATERSHED BOUNDARY DETERMINATION 


The Landsat digital data, as it is received from the EROS 
Data Center, contains a certain amount of spatial distortion when 
compared to commonly used planimetric map projections. Consequent- 
ly, it is necessary to spatially transform either the Landsat data 
to the map projections or vice-versa. For this project the trans- 
formation of the planimetric map data to the Landsat spatial for- 
mat was done for ease of handling. 

Following is a procedure for geometric correction which uses 
watershed boundary determination as an appropriate example of the 
application of this technique. 

The determination of the exact watershed boundary on the 
Landsat digital data is necessary to: 

a. Provide an accurate and meaningful summation of com- 
puter generated data for the watershed; 

b. Minimize computer-related expenses by analyzing only 
the area of the watershed on the Landsat data; and 

c. Provide a useful visual presentation format for the 
graphic display of computer generated data for the 
watershed. 

Procedure 

1. Watershed boundary delineation 

a. Acquire U.S. Geological (USGS) topographic maps in either 
1:24,000 C7.5’ series), 1:62,500 (15* series), or 1:250,000 
scale formats. (Note: The map scale selected should be 
compatible with the digitizing facilities available. In 
the case of this example 1:250,000 scale maps were 

used due to the large size of the watershed.) 

b. Delineate the watershed boundary in pencil on the map(s) 
using the elevation contours to find the ridges and other 
topographic indication of watershed extent such as the 
drainage network (usually shown in blue). Figure 1 shows 
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a small portion of the Middle Fork of the Feather River 
Watershed boundary with the corresponding elevation con- 
tour lines. 

c. Emphasize the final watershed boundary location on the 
topographic map(s) by inking them in using a fine draft- 
ing pen. 

Control point selection 

a. Acquire a 9x9 inch (or larger) simulated color infrared 
Landsat image. Identical to the Landsat digital data to 
which the watershed boundary will be fitted.. 

b. Using the USGS topographic maps (on which the watershed 
boundary has been located) identify a network of points^ 
both inside and outside of the watershed boundaries, 
which can be seen and accurately located on both the 
topographic maps and the Landsat simulated color infrared 
image. Bends in rivers, and lakes and reservoirs will 
probably provide the majority of useful control points 
for most areas. In selecting control points on lakes and 
reservoirs, use the dam or outflow site for point location 
to minimize the effects of horizontal point placement 
error due to the raising or lowering of the bodies of 
water. 

c. Once a complete set of control points has been selected 
and plotted on the USGS topographic map(s), each point 
should be labeled for future reference. Figure 2 shows 
the final control point and watershed boundary location 
for a test watershed, the Middle Fork of the Feather 
River in the Sierra Nevada Mountains of California. 

d. The next step involves the location of the control points 
on the Landsat digital data. The easiest way to accomplish 
this task requires displaying the Landsat data on an inter- 
active display system utilizing a color display and key- 
board. The general area in which the control is located 

is first found by overlaying a Landsat quadrant gxid 
CFigure 3) over a Landsat 9x9 inch image (either in print 
or transparency form) . By using these coordinates to determine 
the general control point location, an image is raised on 
the interactive display system in which the control point 
is located. Using a cursor built into the system or by 
using a grid displayed along with the image on the tele- 
vision screen the operator can find the pixel nearest to 
the control point and record its absolute x and y coordi- 
nates. 









Digitizing the control points and watershed boundary 

a. The USGS topographic map(s) with the watershed boundary and 
control points delineated on it should now be digitized. 

The process of digitizing calculates an x and y coordinate 
value for each control point and for each point on the line 
defining the watershed boundary. Since the coordinates of 
each control point and the line defining the watershed 
boundary are done on the same map base, the actual geometric 
location of the watershed boundary in relation to the control 
points can be calculated and the watershed boundary can be 
located on Landsat digital data. 

b. To complete the digitizing process, a computer compatible 
tape (CCT) of the digitized control points and watershed 
boundary is generated. A suggested format would have the 
tape organized so that the first file would contain the 
control pointls x and y coordinate values, while the second 
Cand last) file would contain the x and y coordinates of 
points on the watershed boundary. 

c. Finally, a list of control points should be prepared show- 
ing Landsat control coordinates (in pixels) versus digitized 
control point coordinates. Table 1 gives an illustration 

of such a listing. 

Transformation of the digitized watershed boundary into Landsat 

coordinates 

a. The first step in the process of defining the mathematical 
relationships between the control points in digitizing and 
Landsat coordinates is to eliminate those control points 
which have been miscalculated in either their digitizing 
coordinates or Landsat coordinates or both. To do this, 
a linear (first order) regression is used in the form of: 

XLandsat = a + b^ Xscanner + b^ Yscanner 

where: Yscanner * 0 

YLandsat = a + b^^ Xscanner + b2 Yscanner 

where; b^^ Xscanner — ♦ 0 

The entire set of control points is processed throijgh the 
first order regression. The control points exhib ■ .ag the 
largest residuals are removed from the data set. -ihe re- 
maining control points are again run through the regression 





Figure 3. Lands at quadrant grid. 
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Digitized Control Landsat Control 


Coorc 

inates 

Coordinates 

Number 

Control Point Name 

X 

y 

X 

y 

1 

Philbrook Reservoir 

0318 

1027 

0867 

1585 

2 

Stirling City Res. 

0096 

1522 

0857 

1761 

3 

Magalia Reservoir 

0001 

1822 

0817 

1898 

4 

Concow Dam 

0005 

.2108 

0925 

1956 

5 

Bardees Bar 

0228 

2120 

1019 

1935 

6 

Cannon Reservoir 

0004 

2846 

0990 

2163 

7 

Oroville North 

0003 

3146 

1009 

2256 

8 

Evans Bay 

0582 

2693 

1316 

2164 

9 

Sky High Bend 

0857 

2398 

1336 

1951 


Table 1 

Listing of geometric control points with digitized coordinates 
and Landsat coordinates. 


» 


equation and the entire process continued until no control 
point has a residual exceeding a preset amount 2.000, 
for example) . Note that since we are fitting points to the 
Landsat format, all residuals are expressed in terms of 
Landsat picture elements (pixels). Care shoald be taken 
during the process of eliminating those control points 
that have large resJ.iduals, so that an "adequate” number 
of points will be left in all areas of the watershed. In 
severe cases, new control points may have to be chosen if 
large areas of the watershed boundary are left unrepresented. 

b. Once the final set of control points has been selected 

(Step 1), the points are processed once more by a quadratic 
(second order) regression equation in the form of: 

2 2 

XLandsat = a + b, Xscanner + b« Xscanner + b„ Yscanner + b, Yscanner 
1 2 3 4 

2 

where: b^ Yscanner + b^ Yscanner » 0 

2 2 

YLandsat = a + b^^ Xscanner + b 2 Xscanner + b^ Yscanner + b^ Yscanner 

2 

where: b^ Xscanner + b 2 Xscanner > 0 

The coefficient values generated by the second order regression 
are now used to transform the watershed boundary (previously 
defined in terms of digitizing coordinates) into a Landsat 
compatible form. 

5. The masking of Landsat digital data with the watershed boundary 

a. One of the more important uses of defining the watershed 
boundary is to "mask” the Landsat digital data. In order 
to do this the computer compatible tape generated in Step 
3.b. of this section is used in conjunction with the 
second order regression coefficients generated in Step 4.b. 
previously to produce a new computer compatible tape in 
which the digitizing coordinate values for the watershed 
boundary and control points have been converted into Landsat 
coordinate values. 

b. This new computer compatible tape of the watershed boundary 
and control points now allows Landsat data in any form 
(raw, recalibrated, classified, etc.) to be masked using the 
watershed boundary. 

Note: The basic procedure described in the fourteen steps above 

represents the geometric correction and subsequent application 
of a specific t^e of ancillary data (planimetric watershed 
boundary) to Landsat digital data. The same basic procedure 
applies to all other ancillary data that is to be "fitted” to 
Landsat satellite data, modified as necessary. 
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APPENDIX II 


REMOTE SENSING-AIDED PROCEDURE FOR 
VEGETATION/TERRAIN ANALYSIS 

INTRODUCTION 

Knowledge of the spatial distribution of vegetation/ terrain 
features within a watershed is necessary for determining evapotrans- 
piration ratios and generating of albedo indices for radiation models. 
Albedo is used directly as an input to the shortwave radiation model, 
which ultimately is used in the watershed-wide calculation of evapo- 
transpiration . 

The procedure used to develop the vegetation/terrain distri- 
bution for a given watershed is based upon the computer-aided 
classification of Landsat digital data. Basically this classification 
involves the selection of training statistics based on the unsuper- 
vised classification of specific test sites located throughout the 
watershed and surrounding area. Those training statistics, represent- 
ing spectral information about specific vegetation/ terrain classes, 
are then used to derive a supervised, watershed-wide classification 
to yield the final vegetation/ terrain distribution. 

Procedure 

1. Gathering of remote sensing and ancillary data 

a. Determine a specific Landsat scene which satisfies 
certain minimum criteria such as : 

1) High sun angle so that shadows on the landscape 
are held to a minimum. 

2) Cloud free imagery in the area to be classified. 

3) Freedom from bad data on the Landsat data such as 
bad scan lines, etc. 

4) Snow-free imagery, to as large an extent as possible 

These criteria imply, then, a late spring through early 
fall date of Landsat imagery. In attempting to select 
a date within this period of time, actual spectral 
response of major vegetation types should be considered. 
If, for instance, annual grassland is a major vegetation 
component within a watershed, then it may be appropriate 


to choose a date in which the vegetation is most highly 
contrasted with the surrounding vegetation complex. 

This date, then, may be in early summer to early fall 
when the grasslands have lost their near-infrared spec- 
tral response due to drying. 

b. Once a date of Landsat imagery has been selected and 
acquired (in both a hardcopy and digital tape form), 
a set of high altitude aerial photographs should be 
acquired that covers the area to be analyzed.. The 
date of this imagery should depict the vegetation com- 
plex in the same physiologic state as the Landsat scene 
selected. It does not necessarily need to be taken in 
the same year as the Landsat scene. 

c. Watershed boundary descriptions such as those used in 
the Geometric Correction and Watershed Boundary Deter- 
mination section should be obtained. These will aid in 
helping to select the appropriate spectral training sites 

2. Selection of spectral training sites 

a. Using the Landsat imagery, the aerial photography, and 
the watershed boundary maps, training sites should be 
selected throughout the watershed and surrounding areas. 
These sites should e^ncompass the spectral variation 
that is likely to occur throughout the watershed. This 
task requires the extensive knowledge of an analyst 
thoroughly familiar with the watershed. 

b. The selected training areas should be located on the 
Landsat digital data and the aerial photography. The 
location of these sites on the Landsat digital data is 
most easily accomplished through an interactive tele- 
vision display of the Landsat data. The x and y coordi- 
nates of the training site block corners can be found 
easily by displaying the Landsat image and determining 
the training site location from its general position 

on either the aerial photography or watershed boundary 
maps. The exact position of the training site should 
be defined by the Landsat image of the training site. 

This area can now be plotted with greater flexibility 
on the aerial photography. This is due to the display 
restrictions usually associated with the Landsat digital 
data. Figure 1 shows the spatial distribution of 
Landsat training sites throughout a sample watershed, 
the Middle Fork of the Feather River. Figures 2 and 3 
illustrate the appearance of a specific training site 
(Block No. 1) on high altitude aerial photography and 
Landsat imagery respectively. 
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3. Unsupervised classification of the training sites. 

a. Once the position of the training sites on the Landsat 
digital data has been determined, a computer compatible 
tape is made with each training site occupying one 
file on the tape. 

b. The next step is to use an unsupervised classification 
algorithm to "cluster" the Landsat pixels at each 
training site into like groups. A commonly used computer 
program for this purpose is the JC^terative Self -Organizing 
Data Analysis T^echnique(A) or ISODATA, which was origi- 
nally developed by Hall and Ball (1965 and 1967), A 
variation on this program, called ISOCLAS, was developed 
by Kan and Holley (1972) specifically for multispectral 
scanner data such as Landsat. This program, ISOCLAS, 

is probably the most commonly used program in use for 
unsupervised classification. Thorough familiarity with 
this, or any other program used for unsupervised classi- 
fication is suggested before attempting this phase of 
the watershed-wide vegetation/terrain classification. 

c. The end product of an unsupervised image classification 
program such as ISOCLAS is a list of image "classes" 
defined by their spectral mean, standard deviation, and 
co-variance. Table 1 shows these parameters for a 
single training site which will be used for the remainder 
of the examples. These classes will now become the basis 
for the final supervised classification for the entire 
watershed. 

d. In order to attempt to assign specific vegetation/terrain 
classes to each cluster, a technique called 7/5 ratioing 
is used to help "spectrally" organize the data. Table 2 
depects the 7/5 ratio information for the training site 
cluster data shown in Table 1. 

e. To help better organize the data for display purposes, 
unique colors are given to each cluster starting with the 
lavenders and the reds for the high 7/5 ratio clusters, 
down through the browns and oranges, the yellows, the 
greens, and finally the blues for the lowest 7/5 ratio 
clusters. A visual output, such as could be seen on a 
television monitor is shown in Figure 4. 

f. Using the aerial photography of the specific training 
site being classified (Figure 1) and the 7/5 ratio, the 
analyst can begin to group and separate the various 
clusters into meaningful vegetation/ terrain classes. 

Table 3 shows the result of this process. Note the arrange- 
ment of the resultant class/cluster groupings. 
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Figure 1. Distribution of Landsat test sites throughout the Middle Fork 
of the Feather River Watershed and surrounding area. 

Area 1 Sierra Valley Mountains 

Area 2 Sierra Valley Area 1 

Area 3 Sierra Valley Area 2 

Area 4 Frenchman Reservoir Area 

Area 5 Quincy /American Valley 

Area 6 Serpentine Area 

Area 7 Davis Lake Area 

Area 8 Gold Lake Area 

Area 9 Silver Lake Area 

Area 10 Middle Fork /Lake Orovi lie 
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Figure 3. 


Landsat "raw" image of ISOCLAS block 1. 
infrared rendition using MSS bands 4, 5 , 


Simulated color 
and 7. 
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Cluster 

MSS Band 4 

MSS Band 5 

!;MSS Band 6 

MSS Band 7 


Mean 

STD 

Mean 

STD 

i 

Mean 

STD 

Mean 

STD 

1 

21 

.27 

1.20 

28.81 

2.26 

33.79 

2.19 

17.34 

1.44 

2 

30 

.37 

1.25 

31.68 

1.75 

35.95 

1,62 

18.19 

1.13 

3 

24 

.27 

2.16 

23.90 

1,63 

29.89 

1,81 

15.99 

1.57 

4 

23 

.73 

1,35 

35.57 

2.08 

38.70 

2,64 

19.62 

2.49 

5 

18 

.97 

2.05 

15.66 

2.17 

31.83 

2.05 

19.21 

1.47 

6 

28 

.34 

1.05 

27.98 

1,65 

31.99 

1.64 

16.44 

1.29 

7 

17 

.53 

1.75 

14.49 

1.72 

26.28 

1.95 

15.42 

1.52 

8 

32 

.29 

1.26 

35.19 

1.86 

39.28 

1.48 

19.77 

1.18 

9 

22 

.33 

2.17 

20.98 

2.33 

36.80 

2.61 

21.57 

1.77 

10 

21 

.79 

1.68 

19.34 

1.74 

27.33 

2.00 

15,17 

1.68 

11 

36 

.61 

5.78 

41.69 

5.81 

43.66 

4,19 

21.67 

2.35 

12 

26 

.99 

2.26 

27.06 

2.15 

39.09 

2,29 

21.82 

1.79 

13 

37 

.52 

1.51 

43.76 

1.88 

47.78 

1.49 

23.98 

1.20 

14 

39 

.94 

3.99 

50.24 

2.93 

54.98 

3.75 

27.24 

1.99 

15 

27 

.13 

4.10 

46.83 

8.59 

46.85 

i 7.28 

24.80 

5.27 

16 

33 

.49 

1.30 

37.30 

2.18 

43.16 

1 1.51 

22.03 

1.52 

17 

26 

.57 

1.72 

38.72 

3.05 

43.94 

!i.98 

i 

22.75 

1.76 

18 

19 

.71 

36.89 

35.86 

35.33 

0.00 

0.00 

0.00 

0.00 

19 


.02 

.15 

4.49 

15.00 

0.00 

0.00 

2.84 

11.05 


Table 1. Cluster means and standard deviations for each of the nineteen 

clusters by Landsat Multispectral Scanner (MSS) bands for ISOCLAS 
block 2. Note the spurious means and standard deviations for 
clusters 18 and 19, indicating that these two clusters will not 
be reliable for further use. 
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Table 2. Band 7/band 5 ratios and color assignments for ISOCLAS block 1 


Band 7/Band 5 Ratio* Numerical Order Cluster Color Assignment 


2.4534 

1 



2.1284 

2 

7 

>' Lavenders § Reds 

2.0276 

3 

9 


1.6127 

4 

12- 


— 

1.5688 

5 

10^ 

• 

1.3381 

6 

3| 

y Browns 8 Oranges 

1.2037 

7 



1.1812 

8 

16j 


1.1751 

9 

ms 


1.1751 

10 

H 

► Yellows 

1.1484 

11 

■1 


1.1236 

12 

ME 


1.1032 

13 

HI 

^ Greens 

1.0960 

14 

K 


1.0844 

15 

HS 


1.0592 

16 


Blues 

1.0396 

17 

HI 



18 

-] 

Black 



19 




*The actual equation for the Band 7/Band 5 ratio is: 

2 X (Band 7 mean scene brightness value of cluster n) 4 (Band 5 mean of scene 
brightness value of cluster n) where n goes from 1 to 17 inclusively. 
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Figure 4. Band 7/band 5 ratio display of ISOCLAS block 1 
based on Landsat data 
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g. Once all of the training sites have been processed in 
this manner, each site should be ground checked to 
verify the classification. Using the Landsat imagery 
thit has been clustered, the aerial photography, and 
topographic maps, each preliminary classification can 
be checked and verified or classification amended, as 
appropriate. Figure 5 shows a display of a final 
unsupervised classification of a single training site. 

4. Selection of training class statistics 

a. Each cluster has a single mean and standard deviation 
in terms of Landsat spectral values assigned to it, 

for each training site. Since it is possible that certain 
clusters could be "confused" with each other, a statis- 
tical test is necessary to determine their ’Separability". 

A methodology for the determination of class separability 
based on the Scheff^ method of multiple comparison has 
been applied to this analysis. Statistical analysis 
using the Scheffe method allows significance probabili- 
ties to be calculated for all cluster pairs. These 
probabilities are displayed in matrix form. Class 
pairs with significance probabilities below specific 
threshold values are considered to be separable for 
classification purposes. Thus individual clusters 
can be statistically compared to other clusters within 
a vege tat ion/ terrain group and all other clusters in 
the vegetation/ terrain classes. Table 4 represents an 
example of the Scheff^ multiple comparison test run 
on all clusters determined to fall into a specific 
vegetation/ terrain class. Cluster pairs determined to 
be separable are designated by the symbol "0". Cluster 
pairs determined to not be separable are designated by 
the symbol "+". The significance level for the parti- 
cular application of the Scheffe multiple comparison 
test was set at a = .10. Clusters form various vegetation/ 
terrain groups (i.e., those failing in some broad band 
of classification; for example, the coniferous forest 
vegetation/ terrain group may include red fir, mixed 
conifer, eastside pine, pure fir, etc.) are systemati- 
cally compared to one another until (if possible) unique 
sets of clusters are found which represent all desired 
vegetation/terrain classes and are all separable using 
the Scheffe test, 

b. The resultant set of clusters now represents each 
vegetation/ terrain class desired for the watershed- 
wide classification. Table 5 lists the mean, standard 
deviation and co-variance for each selected cluster. 
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Figure 5. Classified ISOCLAS display of block 1 based on 
Lands at data 
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Table 4. Scheffe multiple comparison test for brush and brush- 
related vegetation types. "0" indicates no confusion 
while indicates confusion between class pairs at 
a significance level of .10. 
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Cluster 

MSS Band 4 

MSS Band 5 

MSS Band 6 

MSS Band 7 

Mean 

STD 

Mean 

STD 

Mean 

STD 

Mean 

STD 

1 

30.37 

1.25 

31.68 

1.75 

35.95 

1.62 

18.19 

1.13 

2 

28 .34 

1.05 

27.98 

1.65 

31 .99 

1.64 

16.44 

1.29 

3 

21 . 79 

1.68 

19.34 

1.74 

27.33 

2.00 

15.17 

1.68 

4 

36.61 

5.78 

41.69 

5.81 

43.66 

4.19 

21.67 

2.35 

5 

34.67 

1.55 

39.02 

1.55 

43.65 

1.32 

22.81 

1.19 

6 

21.10 

2.14 

14.45 

1.41 

64.66 

1.82 

42.94 

1.68 

7 

34.34 

1.53 

39.90 

1.86 

48.43 

1.64 

26,13 

1.36 

8 

36.45 

1.54 

42.43 

1.56 

46.60 

1.99 

24.43 

1.66 

9 

24.95 

2.29 

25.83 

1.95 

44.22 

2.23 

27.54 

1.57 

10 

31.60 

1.39 

34.68 

1.92 

38,97 

1.79 

20.46 

1.66 

11 

21.73 

2.12 

20.13 

2.44 

20.10 

1.39 

25,81 

1.04 

12 

22.70 

2.56 

21,86 

2.43 

23.13 

2.25 

10.80 

1.83 

13 

17.13 

1.66 

11.88 

1.80 

28.96 

1.37 

17.63 

1.13 

14 

17.36 

1.71 

12.43 

1,72 

32.91 

1.52 

20.64 

1.07 

15 

26.49 

2.28 

26.49 

2.29 

43.90 

1.61 

26.06 

1.77 

16 

17.98 

1.98 

12.48 

1.23' 

43.69 

2.53 

28,63 

2.35 

17 

17.47 

1.72 

13.32 

2.04 

31.41 

1.51 

19.68 

1.35 

18 

18.81 

2.08 

13.63 

1.90 

51,63 

3.17 

34.30 

2.27 

19 

32.23 

3.60 

38.82 

2.81 

43.63 

2.91 

22.92 

2.38 

20 

15. '72 

1.69 

10.00 

.84 

6.25 

1.23 

1.07 

.65 

21 

16.30 

.1.74 

11.92 

1.71 

23.71 

1.76 

14.37 

1.41 

22 

14 . 73 

1.17 

6.73 

.49 

4 . 79 

.93 

1.04 

.66 

23 

17.78 

1,78 

13.31 

1.73 

30.67 

1,32 

19.21 

1.16 

24 

19.39 

1.91 

14.63 

1.79 

45.67 

1.74 

29 . 84 

1.60 

25 

16.09 

1.75 

9.58 

2.09 

11.15 

2.12 

4.91 

1.98 

26 

14.16 

1.84 

7.73 

6.75 

4,97 

1.17 

1.19 

,67 

27 

22.86 

1.78 

19.73 

1.68 

37.34 

2.10 

21.64 

1.38 

28 

19.39 

2 . 12 

14.22 

2.14 

31.76 

1 . 24 

18,86 

1.11 

29 

25.53 

2.30 

24.49 

2.06 

40.65 

1.77 

23.01 

1.32 

30 

28.89 

4.24 

25.85 

4,44 

29.48 

3.60 

13.77 

2.83 


.Table S. Cluster means and standard deviations for each of the thirty 
clusters by Landsat Multispectral Scanner (MSS) bands for 
the final CALSCAN training set. 
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Supervised classification of vegetation for the entire watershed 

a. The cluster means, standard deviation, and co-variances 
for the set determined above are now used to drive the 
final, supervised classification of the watershed using 
digital Landsat spectral data. The computer compatible 
tape of the 'watershed boundary (mask) generated in 
Step 14 of the Geometric Correction and Watershed 
Boundary Determinations section (Appendix I) is now 
used to control the area on the Landsat data that is to 
be classified. The ’’raw" Landsat digital data in con- 
juction with the "mask” tape are now run through a 
supervised classification routine. The classification 
routine commonly used in conjunction with the classifi- 
cation of Landsat data is a discriminant analysis routine 
called LARSYSA. LARSYSA is a classification algorithm 
based on a maximum likelihood classifier developed at 
Purdue University. A similar version, based on LARSYSA, 
used at the University of California, is called CALSCAN. 

Both classifiers perform the same basic function, assign- 
ment of each pixel within the watershed to one of the 
training clusters selected from the unsvipervised classi- 
fication of the training sites. The final result for 
such a classification can be sef;n for the ISOCLAS test 
block in Figure 6, while the final watershed-wide classi- 
fication is shown in Figure 7. 





Figure 6. CALSCAN discriminant analysis classfication of block 1. 

The statistical training for this classification was 
derived form the ISOCLAS classification and applied to 
the Landsat raw digital data. 
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Figure 7. CALSCAN discriminant analysis vegetation/terrain 
classification of the Middle Fork of the Feather 
River Watershed. 
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APPENDIX III 

STEPWISE PROCEDURE FOR TOPOGRAPHIC ANALYSIS 


Topogra^phic characteristics of the watershed are some of the 
major factors affecting the amount of solar radiation. These 
topographic characteristics include elevation, slope, and aspect. 
The following is a procedure for quantification of these parameters 
along with three different methods of control point selection and a 
set of algorithms for elevation interpolation and slope and aspect 
generation. 

For a large portion of the country, digital terrain computer 
compatible tapes are available. These tapes may be obtained 
either from the Defense Mapping Agency (DMA), Topographic Center, 
or from the U.S. Geological Survey (USGS) , Digital Elevation Model. 
If these digital tapes for the study area are available, the part 
of the following procedure for digitizing the elevation from USGS 
topographic maps may be eliminated. Otherwise, the following pro- 
cedure will be adequate for preparation of elevation, slope, and 
aspect maps and digital tapes required for resource inventory 
and management . 

1. Procedure for topographic analysis 

a. Acquire USGS topographic series maps, scale 1:250,000 
so that the area analyzed is covered. (Note: the same 
maps used for the Geometric Control and Watershed 
Boundary Determination Section may be used for this 
phase of the work as well.) 

b. If the topographic maps used do not have the study 
area outlined on them, outline the area at this time. 

c. If there are two or more topographic maps needed to 
cover the study area, carefully trim and join the 
maps, taking care to maintain the planimetric accuracy 
as much as possible. 

i • . ' 

d. If I the topographic maps used do not have control 
points (for geometric correction purposes) located 
on them, see the Geometric Correction and Watershed 
Boundary Determination Section (Appendix I) of the 
Manual. 
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e. The selection of topographic control points is the 
next phase of the topographic analysis procedure. 

The actual technique used depends on the nature of 
the elevation, slope, and aspect generating algorithms 
used. Three methods for selection of control 
points are discussed in Section 2. A series of 
algorithms for elevation interpolation and slope and 
aspect generation is described in Section 3. 

f . Each topographic control point is given x, y and z 
coordinate values. The x and y values are relative 
to the geometric correction points and are used for 
two-dimensional spatial location of the z or eleva- 
tion value. The topographic control point and geo- 
metric control point networks, as located on the 
USGS topographic maps, are now digitized. This 
process assigns x and y values to all control 
points. The geometric control points should be 
digitized in a predetermined order so as to facili- 
tate their accurate identification. The topographic 
control points should be digitized in such a fashion 
as to facilitate the assignment of a specific z 

(or elevation) value, to each control point. 

g. To complete the digitizing process, a computer 
compatible tape (CCT) is generated for the topo- 
graphic control points. The z value for each control 
point should be added to the digitized topographic 
control point data. This may be done on either the 
CCT or on punched cards (the x and y coordinates 
coming from the CCT with z values punched in manually). 

h. Using the regression coefficients generated pre- 
viously or that will be generated specifically for 
this analysis (for both counts see the Geometric 
Correction and Watershed Boundary Determination 
Section, Appendix I) in conjunction with the x, y 
and z values fof each topographic control point, 
the x and y coordinates (in digitizing coordinate 
values) are transformed to Landsat-based coordinates. 

i. The transformed x and y coordinates as well as the 
unchanged z value are now processed through an 
"interpolator*' to assign a z (elevation) value to 
each Landsat pixel position. This "interpolation" 
of intermediate x y coordinate positions for eleva- 
tion can be accomplished by a number of algorithms. 

In Section 3, one type of algorithm for elevation 
interpolation is described. , 


I II -2 




j . The interpolation is processed by a computer, the 
output usually being in the form of a CCT and supple- 
mentary paper output. The elevation interpolation 
may be displayed using the ~)CT on an interactive 
television display system. An example of the final 
product, using different colors for each elevational 
zone, can be seen in the case study. 

k. ITie generation of slope and aspect maps is done by 
using another algorithm which uses the interpolated 
elevation values as the input data. An example of 
the output products of this algorithm can be found 
in the case study. 

2, Selection of topographic control points 

Three procedures will be shown in this section to 
give the user some idea of the possibilities for 
topographic control point selection. The ultimate 
decision for location of topographic control points 
will depend on the available computer system (and 
associated software routines) as well as the nature 
of the algorithms used for generating the elevation, 
slope, and aspect data. 

Method a. This procedure involves the use of elevation 
contour isolines. Each isoline falling within the 
study area is digitized. By digitizing one eleva- 
tion zone at a time, the data may be organized so 
that assignment of elevation values to digitized 
isolines is very sijnple and straightforward. If a 
computer compatible tape is used to store the digitized 
contour isolines, a separate file may be generated 
for each different elevation zone. Figure 1 shows 
a hypothetical distribution of control points 
using this technique. It should be noted that as 
few points as possible should be selected along each 
isoline to keep processing costs to a minimum. This 
method of control point selection is probably the 
most accurate and undoubtedly the most expensive 
due to the large number of control pointsusually 
selected. 

Method b. The selection of topographic control points 
along ridge tops and valley bottoms is a commonly 
used procedure. Control points are located on ridge 
tops (including watershed boundaries, if they are 
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to be used) and in the valley bottoms (generally 
along stream courses) . Control points are selected 
at the ridge intersections and midway (approximately) 
down ridges running into valleys. Subsequently, 
control points are also located at stream inter- 
sections and midway up streams from their head- 
waters (Figure 2) . The control points selected are 
digitized. The difficulties of this method involve 
the location of the ridge tops and valley bottoms, 
and the assignment of elevation values to each 
control point . 

Method c. This method uses a predetermined grid which 
is laid over the study area. Each grid intersection 
is located as a control point and the elevation recorded. 
This system does not necessarily require the use of a 
digitizer due to the systematic nature of the control 
point network (however, a iligitizer can be used) . Again, 
the assignment of an elevat ion value to each control 
point is timely; however the nature of the control 
point selection procedure allows a tight control over 
the number of control points selected, thus helping 
to reduce costs. 

Note: In order to assure as accurate values as possible for elevation 

slope, and aspect, control point selection should occur not only within 
the study area but also just beyond the boundaries. This assures 
that values along the actual study area boundary are reliable. 

3. Sample algorithms used for elevation, slope, and aspect 
generation 

The following formulas are those used by the Remote 
Sensing Research Program at the University of California 
at Berkeley, for topographic analysis. They are 
representative of those commonly in use. 
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Figure 2. Location of control points on ridge tops and valley bottoms, 
Note the location of occasional control points along long 
stretches of ridge tops and valley bottoms. 
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The digitized values of control points are used to drive the 
elevation-determining algorithm, 

where j Ai = the distance from (x,y) to control point i, 

I " I 

Ei = the elevation of control point i, 

n = the number of control points within range 
of (x,y) 

E(x,y) = the elevation of point (x,y) . 




which uses the x,y,z coordinates of the control points to assign a 
z coordinate (elevation) to each Landsat pixel. The actual z 
coordinate represents a zone of elevation per pixel rather than an 
absolute value, 

r"“'" 

Determination of a "normal" vector is necessary for pixel 
slope and aspect calculation. This "normal" vector is a directional 
variable that originates in the center of each pixel and locates a 
direction that is 90 degrees to the surface of the pixel. This 
vector is expressed in triple coordinates (n^, ny, rij,) , where the 
center of the pixel is considered to be the origin (o, o, o) , and 
the triple coordinates are the coordinates in three dimensions 
(x, y, z) . 

For determining slope and aspect, the "normal vector for a 
pixel is computed by using the pixel's position relative to the 
pixels to the right and below. The "normal" vector, ii, is then the 
cross product of the two vectors in x (the pixel to the right) and 
y (the pixel below) directions. 


pixel to the right 




Using the vector’s cross-product formula, the normal vector is 

n = (x + 1, y^ z ^ z') x (Xi y + 1, z"- z) 

=». |[(x + l)(y+l) - xy] - [(x+l)(x"-z) - x (z - z]» 

[(n + l)(y+l) - n y]j 
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Once the ’’normal" vector, n, is computed, the slope is calculable 
since it is the angle between the normal vector and the vertical 
vector Co, o, 1), 0, subtracted from jin radians. Therefore 

G = arc cos { ) and 

lUll 

7T 

slope " y “ 9 ill radians, where 

II Sll = 

then slope is converted to percent. 

The aspect is determined from the (x,y) projection of the normal vector. 
The degrees away from true North will determine the aspect, 

n 

<}) = arc tan (~^) , where 
n 

X 

<f) is the angle from the x-axls in radians. 

This is then converted to degrees and added to 90 degrees to give the aspect 
in degrees from the true North direction. 


I 
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APPENDIX IV 

PREPARATION OF CLIMATIC DATA MAPS FOR USE IN 
RADIATION AND EVAPOTRANSPI RATION MODELS 

This appendix entails three sections: 1) collection and inter- 
polation of climatic data based on ground meteorological stations 
data; 2) acquisition and processing of NOAA-4 satellite thermal data; 
and 3) a complete version of a computer algorithm for geometric cor- 
rection and temperature calibration for NOAA-4 Satellite data. Sec- 
tion 3 was originally provided by NOAA’s National Environmental 
Satellite Service and was modified by UC Dayis, Department of Elec- 
trical Engineering (Algazi group). One may develop his own algorithm 
based on the technique described in Section 3 and geometric package 
programs available from the NOAA-NESS office. 

1. Collection and interpolation of data from ground 
meteorological stations 

a. Identification and selection of ground meteorological 
stations 

Subsequent to the delineation of the watershed 
boundary (described in Appendix T) on a useable map 
base (preferably a USGS quadrangle, the scale of 
which will depend on the size of the watershed) one 
should consult the publication titled "Climatological 
Data" provided by the National Oceanic and Atmospheric 
Administration's (NOAA) Environmental Data Service. 

This publication is provided on a monthly basis, 
state-by-state, and includes daily and long term 
average meteorological information with varying 
degrees of specificity. Provided within the pub- 
lication is a map of the particular state that is 
of concern. All of the stations that provide data 
to the Environmental Data Service are indicated, as 
are the types of information they collect. 

With this information, the number and type of 
stations within the study area may be selected. 

b. Registration of the ground meteorological station data 

In order to prepare the climatological informa- 
tion for use in a Landsat based study, the information 
must be made compatible to the satellite data coordi- 
nate system. To accomplish this, the meteorological 
stations being used should be plotted on an acceptable 
map base (preferably a USGS quadrangle). Next, the 
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Universal Transverse Mercator (UTM) coordinates for 
each station should be determined from the map. Now, 
through usd of the regression equations developed 
during the geometric correction of the original ' 

Landsat data, the Landsat coordinates may be calcu- 
lated for each ground meteorological station. The 
climatological data is now registered to the satellite 
data base and may be used in direct conjunction with 
raw Landsat data or other data sources based on the 
Landsat coordinate system. 

c. Recording and coding of ground meteorological station 
data 

The climatological data derived from each of the 
ground meteorological stations, as well as important 
information such as date of data acquisition, station 
name, UTH coordinates, Landsat coordinates, and eleva- 
tion, should be noted and recorded in a useable format. 

A convenient way to accomplish this is by the use of 
computer cards. Under a format developed by the user, 
all information collected about each station may be 
recorded on one or more computer cards . The information 
is now in an acceptable format for use in a data base 
computer routine such as RSRP’s MAP IT program or an 
interpolation algorithm, 

d. Climatic data interpolation 

Since climatic data representing the entire water- 
shed is collected from a number of point sources, there 
is a need to expand this data base so that climatic 
information will exist for each Landsat pixel (picture 
element). To accomplish this, an interpolation algorithm* 
was developed at RSRP that will mathematically generate 
data values for each pixel based on the ground station 
data (i.e. temperature, humidity, precipitation, 
evaporation, etc.). Upon completion of interpolations 
for each of the climatic variables of concern, the 
data is ready to be used as input into data base 
programs such as RSRP's MAP IT, or simulation models 
for short-wave, long-wave, or net radiation and evapo- 
transpiration estimation. 


*See following page for interpolation algorithm equation. 


*The interpolation algorithm is described by the following equation: 


s — 
i=l(D.)^ 
D(x,y) = — - — — 
n 1 

^ ^2 

1=1(0^) 


where : 

D. = the distance from the Landsat pixel 
^ being considered (x,y) to the ground 
meteorological station (i) expressed 
in terms of Landsat pixels. 


D. = the data value associated with the 
^ ground station. 


n = the number of ground station locations 
within range of the pixel being con- 
sidered (x,y) , 


D(x,y) = the calculated data value of the pixel 
Cx,y) based on n ground station 
locations . 


The derived value D(x,y) represents the two coordinates of a three 
dimensional (x,y,z) computer database (such as MAP IT) , registered 
perfectly to the x,y Landsat coordinate system. 


e. Use of pz'epared climatic data, in simulation models 

After registration and interpolation of the climatic 
data, the infoimation is ready for use as input into a 
number of simulation routines. Among these are models 
for estimating short-wave, long -wave, and net radiation 
and a number of models for evapotranspiration estimation. 
The output of such models is written onto magnetic 
tape and line printer class maps . The digital data on 
magnetic tape may be used to generate a visual pro- 
duct such as with RSRP’s film writing device known as 
IGOR (Image Ganged Optical Reproducer), or as input 
into higher level simulation models, such as the water 
yield model developed by the joint Federal-State River 
Forecast Center. 

f. Generation pf final products 

Upon completion of the data preparation and model 
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simulation activities, there is often the need to 
transform the derived data into a presentable form. 
Digital information, as it is derived from the solar 
radiation and evapotranspiration models, may be used 
as input into a film writing device such as IGOR. 
Selection of a logical sequence of colors to repre- 
sent the different classes of digital information 
aids in the production of a meaningful and visually 
pleasing product. Following this, concise descrip- 
tive labels and a statistical summary should accom- 
pany each visual product to aid in its interpretation 
and effectiveness. 

NOAA-4 satellite thermal data processing 

Under the present climatic stations network, areal 
distribution of temperature for radiation models may be 
obtained adequately from ground meteorological stations. 
Data obtained form the NOAA-4 satellite offers an al- 
ternate method for generating a temperature map for a 
large area. The following text gives a) a brief summary 
of relevant information about the NOAA-4 satellite] b) 
an outline of required computer system facilities for 
implementing this procedure] and c) a step-by-step pro- 
cedure developed by the Remote Sensing Research Program 
(RSRP) for NOAA-4 thermal data processing utilizing 
digital tapes obtained from the National Environmental 
Satellite Service (NESS) of NOAA. The algorithm used 
for part of the data processing (geometric correction 
and temperature calibration) has been provided by NESS 
and modified by grant participants in the Department 
of Electrical Engineering at UC Davis. 

a. A brief description of the NOAA-4 satellite 

The NOAA-4 is an environmental satellite in a 
sun- synchronous near -polar orbit. This satellite 
is an evolution of the NOAA series satellites launched 
by NASA and administered by NOAA after becoming 
operational. It is equipped with various sensors 
and scanners for obtaining meteorological data. 

The data used in this procedure was obtained by 
a Very I^gh Resolution Radiometer (VHRR) , a two- 
channel scanner which detects energy in the visible 
spectrum (.6 to .7 micrometers), as well as the 
thermal infrared (10.5 to 12.5 micrometers). The 
satellite orbit has an inclination of about 102°, 


an altitude of about 14S0 km, and a period of 115.0 
minutes . Its daytime pass over the earth is in a 
general Northeast to Southwest direction crossing 
the equator at an angle of 78® at approximately 
9:00 AM local solar time . The VHRR scans west to 
east, from horizon to horizon,, at the rate of 400 
lines per minute. Its spatial resolution at the 
satellite subpoint is about 0.8 km. 

b. Required system capabilities 

The following computer facilities are required fo.r 
implementation of the procedure: 

1) A magnetic tape unit 

2) A disk pack 

3) A video display system capable of displaying 
3-bit data as a black-and-white image, with 
some means of obtaining the coordinates of 
points in the image. 

The software requirements include: 

1) A program for reading the NOAA tape and 
generating a disk file for corrected or un- 
corrected imagery (See Appendix V) 

2) A program for displaying the image on a tele- 
vision monitor 

3) A means for obtaining the relative coordinates 
of a point in the displayed image 

4) A set of algorithms for implementation of 
various required; m.athematical equations on 
the data file 

5) A program for linear regression performance 

6) A program to produce a histogram of the data 
values 

7) A program to transform the data files to the 
desired coordinates 

With the exception of a), which is given in Appendix V, 
almost all of the software programs should exist witliin 
any software package for an image analysis and display system 


c. NOAA data acquisition 

The computer compatible tapes for both visible 
and thermal data can be obtained from NESS headquarters 
in Washington, D.C. The data for this study was ob- 
tained form grant participants at UC Davis who had 
copied 1600 bpi NOAA tapes to 800 bpi tapes, with a 
slight change in format. A written description of 
the important information for geometric correction 
was received along with the tapes . This information 
included satellite altitude, midpoint of the image, 
format of the tape, and the raw coordinates of one or 
two reference points. 

d. Location of the area of interest on the raw thermal 
data : 

1) If possible, estimate coordinates of the area of 
interest from the coordinates of a known reference 
point given in the written description. 

2) Create a raw (uncorrected) file which includes 
the desired area. For the first trial, it is 
suggested that the image be reduced by a factor 
of four or more, and the starting coordinates 
adjusted accordingly, to allow a large margin 
for error. 

3) Display the file and locate the area of interest. 

If it is not in the display, re-estimate start- 
ing coordinates and repeat Step 2.b. 

4) Determine the coordinates (pixel and line numbers 
with respect to the beginning of the tape) of 
the upper left corner of the area of interest , 
and also determine the size of the area. 

5) Using the new coordinates, repeat Steps b and c 
without reducing the image. 

6) On the new display, choose a reference point, 
preferably near the center of the image, which 
can be accurately located on a map. Determine 
its coordinates with respect to the tape file, 
and determine its exact latitude from a plani- 
metric map. Estimate the relative coordinates 
of the reference point (pixel) in the desired 
geometrically corrected output image. 

e. Creation of a geometrically corrected image of the 
study area 
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Use the parameters determined in c.2. and the 
sate Hite height and period data to create a geo- 
metrically corrected image. The meanings of the 
parameters required by the NOAA program for geo- 
metric correction are as follows: 

SIZEX - line length of output image in pixels 

DELTAX - sampling distance along a line of input image 

SIZEY - number of lines in output image 

DELTAY - sampling rate of lines in output image 

RAWXl - x-coordinatb of reference pixel with respect 
to beginning of tape file 

RAWYl - y-coordinate of reference pixel with respect 
to beginning of tape file 

HEIGHT - height of satellite in kilometers at time 
of image acquisition 

PERIOD - period of satellite (115.0 min.) 

REFLAT - latitude of reference pixel 

IMAGXl - relative x- coordinate of reference pixel in 
output image 

IMAGYl - relative y-coordiante of reference pixel in 
output image 

A complete description of the geometric correction algorithm 
is given in Section 3 of this appendix. 

Creation of a histogram of the data values on the geo- 
metrically corrected image of the study area 

Inspect the histogram to determine the range of 
valid data values. A temperature range of 0° - 40°C 
would correspond roughly to a data value range of 120-200. 

NOAA data enhancement 

The NOAA satellite data will probably have some 
visible banding (bad data lines primarily due to 
temperature variations of the sensor at the time of 
scanning) and "noise" (extremely high or low values 
inconsistent with all surrounding pixels) . These 


two problems are handled separately. The corrected 
data is in each case generated as a new band appended 
to the data file. The banding is assumed to be the 
result of variation in the IR sensor, which could 
cause data values in adjacent scan lines to be in- 
consistent with each other, i.e. the same values 
could represent different ground temperature values. 

Values within one line are assumed to be consistent, 
since the amount of time required for one scan of the 
3 tudy area is small compared to the time between scans , 
Since each scan line cuts across a diverse terrain, 
the line averages are expected to vary smoothly, and 
sharp deviations are probably the result of sensor 
variations. 

The procedure developed for reducing banding 
effects is the following: 

1) Compute the mean of all pixels in line 

n with Values within the range pf normal data. 

2) For each line (n) , compute an adjusted 
mean A^, obtained from the sequence by 

"smoothing”. For this project, was taken 

to be simply an average of the means of the 
five lines centered on n, weighted by binomial 
coefficients: 

2 

\ ^ S,k+3^n+k 

3) Increment the value of every point in line n 

by the interger closest to • In the 

data to which this procedure was applied, there 
was only one change of 3 units, about 7% of 
the lines were changed by ^2, and about 45% 
by +_1. A change of 2 units is approximately 
equivalent to one degree of centigrade temperature 

There may be some apparent noise in the image. The 
following procedure was used to correct for noise. This 
technique may be applied by using either vertical or 
horizontal lines in the image. For each point in the 
smoothed image, the following is applied: 

Ij Obtain its smoothed I R value (Vj from the data 
file 
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2 ) Obtain the values U and W which are the smoothed 
IR values for the closest preceding and following 
pixels, along a vertical or horizontal line, whose 
values are within the range of valid data. 

3) If V is not a valid data point, or if V is more 
than some tolerance (T) outside the range of 
values determined by U and W, then replace V by 
(U+W)/2. The process does not depend critically 

on the tolerance T, unless it is too small or large. 
T=10 appears to be reasonable. 

Determination of voltage wedge values and the space 
view average 

The voltage wedge values and space view average 
are used for transforming the raw values to voltages. 
Compute average values using the header data for all 
the lines in the study area. The space view average, 

S, is in byte 12, and the wedge values W(0) , . . . ,W(5) 
are in bytes 13-18. Also note the IR value (byte 20) , 

T1 (byte 21), and T2 (byte 22), for later use. Typical 
values of these parameters are: W(0)=216, W(l)=i81,v 
W(2)=144, W(3)=108, W(4)=70, W(5)=35, S=30, IR=150, 
Tl=180, T2=186. See Appendix V for further details. 

Transformation of raw values to voltage values 

Tills is a special procedure developed for data 
acquired prior to September 1975, for which some of 
the calibration data was recorded incorrectly. Com- 
pute the voltages V(0) , . . , ,V(5) corresponding to the 
wedge values W(0) , . . i ,W(5) by V(I) = (W(5) -W(4)) + 

5.8 + (1-5.0). To compute the voltage corresponding 
to a data line value, interpolate between the two 
closest voltage values. 

Transformation of voltage values to temperature 
values 

Following the procedure described in Appendix 
V, compute the values and T^. On the "Radiometer 

Radiance Calibration Curves" graph paper, plot 5.8v 

and O^K, and V ,T . Connect the two points by a 
c c 

straight line. Carefully read off 10-15 points in 
the range 0-40°C. Run a linear regression, taking 
as the dependent variable temperature, and as the 
independent variables voltage V, V , 1/(5.8=V) and 
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ClnC5.8-V)) ■ . This choice of variables gives 
a good fit, with residual values of 0.1°C or less, 
except for an aberration in the graph at about 
12-16°C, where residual values may be 0.2-0.3^C. 

Use the above-mentioned equation to generate tem- 
perature values in °C, 

k. Transformation of temperature values to the desired 
coordinates 

If it is necessary to transform the temperature 
map to different coordinates, the following procedure 
is recommended: 

1) Display the temperature band on the display 
system 

2) Locate likely control points and determine 
their coordinates fairly accurately. This 
can be quite difficult on NOAA data. Most 
good points will probably be on bodies of 
water. 

3) Print out the temperature values for blocks 
(about 15x15] of pixels centered on these 
control points. Comparing the printout 
with a map or a displayed image of the area 
in the desired coordinates, trace out contour 
lines (e.g. lake shorelines) to locate the 
control points more precisely. 

4) Using the pairs of coordinates determined 
above, run a linear regression between the 
NOAA coordinates as dependent variables, and 
the desired coordinates and their squares as 
independent variables, 

5) Transfer all temperature values to the 
desired coordinates using the relationship 
established in Step d. 

3. A computer algorithm for geometric correction 

The following computer algorithm (see Figure lV-1) 
may be used for geometric correction. The temperature 
calibration is done according to the procedure developed 
by NOAA-NESS. A modified version of this procedure, 
developed by the personnel at UC Davis, is described on 
the pages following. 
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a. NOAA thermal calibration 

NOAA has provided us with a procedure for con- 
verting the IR data into surface temperature, but 
for images acquired prior to September 1975, the 
required calibration data was incorrectly recorded 
on the header records, so the procedure does not 
work. The following is a modified procedure which 
has been shown to be reliable by the Department of 
Electrical Engineering at UC Davis . 

b . Procedure 

1) Determine the voltages corresponding to the 
voltage wedge values on the data records (see 
NOAA Tape Format, attached) . Denoting the 
wedge values by Wi (i=0,.,.,5) with corres- 
ponding voltages Vi, and denoting the space 
view average by S, the relationship is: 

Wq-S 

Vi = + 5.8 + (i-5) 

W5-W4 

Prepare the voltage wedge table: 


Voltage Wedge 
Values 

Voltage 

Wo 

Vo 

Wi 

Vl 

W2 

V2 

W3 

V3 

W4 

V4 

W5 

V5 


• 2) Determine the voltages of the thermistors Ti 
T .2 by interpolating between values in the 
voltage wedge table. This will establish 
voltages Vpi and V'i’2. 

3) Determine the temperature of the IR calibra- 
tion target, denoting the temperature indi- 
cated by Ti as Tpi, and by T2 as T-p2* 
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T-ri = 3.60 Vji + 11.40OC 
T.j-2 = 3.57 Vt2 + 11.68°C 


Use the average of these values as the calibra- 
tion temperature T^. 

4) Determine the voltage corresponding to the IR 
calibration target value using the voltage wedge 
table. Denote this value hy V^, 

5) Plot on the graph labeled "Radiometer Radiance 
Calibration Curves" a point at 5.8v and 0°K, and 
a second point at Vq and T^. Connect these two 
points by a straight line. 

6) A data value is now converted to temperature by 
first determining the voltage by interpolating 
between values in the voltage table and then 
finding the temperature from the graph. 

c. Thermal calibration example 


The record header for 5-24-75 is: 


Character: 

12 

13 

14 

15 

Value (hex): 

IE 

D8 

B5 

90 

Value (decimal) : 

30 

216 

181 

144 

Step #1 : 

Vi 

Ws 

-S 

+5.8 

"Wc- 

-w. 


16 

17 

18 

19 

20 

21 

22 

6C 

46 

23 

64 

96 

B4 

BA 

108 

70 

35 

100 

150 

180 

186 


Ci-5) 

0.8 


= i + 0,657 


The voltage wedge table: 


Wedge Values 

Voltages 

Wo=216 

Vo=0.657 

Wi=181 

Vi=1.657 

W2=144 

V2=2.657 

W3=108 

V3=3.657 

W4=70 

V4=4.657 

Ws=35 

V5=5.657 
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step »2: Vti =. + Vj = » 1.657 = 1.684v 

W **T 

''T2 = * Vo = t 0.657 = 1.514V 

Step #3: T-ri = 3.60 V.j.j + 11.40 = 17.46°C 

T^2 = + 11-68 = 17.08OC 

Tc = %(Tti + Tt2) = 17.30C 
Step #4; Vc ^ 1.657 = 2.4rv- 

Step #5: See attached graph (Figure IV-l) 

Step #6; If the radiometric data value (D) is 188: 

„ _ 216-188 . „ , ,, 

“ 216-181 6.657 - 1.46v 

From the graph, Tp = 36.8°C 
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NORA - R ROUTTKE TO GENERRTE THE RRW OR QEOMETRICRi.Lt' 

CORRECTED o;STRETCHED> NORA IMAGE FILE FROM DflTfl 
RFAD FROM A MORfl TAPE. STRETCHING CORRECTS 
FOR PRNORRMTC <EARTH CURVATURE? AND 
SKEW <ERRTh i;OTAl lOtOERRORS. 

RSRP. U. C. BERKELEY 6/1-1/7S 

ORIQINAI..LV ADAPTED FROM A PROGRAM BV U. C. DAVIS, 

ELECTRICAl ENGINEERING DEPARTMENT 

PARAMETERS : 

INPUT 

■ SI7FX = NUMBER OF PIXELS ALONG A LINE IN THE OUTPUT IMAGE FILE 
DE1..IAX = SAMPLING DISTANCE ALONG A LINE OF THE INPUT DATA 
S17EV --- NUMBER OP LINES IN THE OUTPUT IMAGE 
DEt.TAV = SAMPLING RATE OF LINES IN THE INPUT IMAGE 
HFKTHT = HEIGHT OF SATELLITE IN KILOMETERS 
PERIOD - ORBITAL PERIOD OF SATELLITE IN MINUTES 
RAMX;1. RAW REFERENCE PIXEL NUMBER ALONG A LINE OF INPUT IMAGE 
RAMVd RAW REFERENCE PIXEL LINE IN INPUT IMAGE 
REFLHT = LATITUDE OF REFERENCE PIXEL IN DEGREES 
IMAOXd - LOCATION OF RAWXl IN OUTPUT IMAGE 
IMAGVl = LOCATION OF RAWVl IN OUTPUT IMAGE 
RAW - 0 TO GENERAIE RAW IMAGE FILE E< RAWXl, RAWVl? WILL BE 
LOCATED AT <1,1? IN OUTPUT FILE! 

1 TO GENERATE A GEOMETRICALLY CORRECTED IMAGE FILE 

PROGRAM 

RRWX0 = RAW PIXEL HUMBER OF THE FIRST PIXEL IN THE FIRST LINE 
OF THE INPUT IMAGE USED FOR COMPUTING OUTPUT PIXEL 
<1,1?, RAWX0 NEED NOT BE AN INTEGER. 

RAPJV0 --- RAW LINE NUMBER OF THE FIRST LINE USED FOR COMPUTING 
OUTPUT PIXEL <1, 1? 

INTX0 r.- corrected pixel NUMBER OF RAWX0 BEFORE TRANSLATION 
TO OUTPUT IMAGE - 

INTXr). =-■ CORRECTED PIXEL NUMBER OF RAWXl BEFORE TRANSLATION 
TO OUTPUT IMAGE 

SVNCXl? = FIRST SVNC MARK FOUND 

SVNC<2) - SECOND SVHC MARK FOUND 

MID PIXEL NUMBER OF MIDLINE OF INPUT IMAGE 

IHFTA - SATELLITE LATITUDE 

DEI, PHI INCREMENTAL GEOCENTRIC ANGLE 

DEHETA = INCREMENTAL GEOCENTRIC SKEW ANGLE 


TNCl.LIDE "NOAASPEC" 

DATA PI/3. 141S93X 

GET PARAMETERS AND OPEN FILES 

CALL INPAR 


Figure IV-1. Computer algorithm for geometric correction. 
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c 

:?0 


c 

50 

60 

70 

C 

80 

SO 

100 

i.l 0 

130 

C 


CttI I.. irjpEN 

SET UP PPMtDVS HEADER 
DO 30 I=0.- 255 

hf:hiki>=0 

H£HI''<0>-1 
HEHI,H1>--1 
HtHIKJ0>?=l 
HERD<li)=l 
HEHIK12>=SIZEX 
HEHIX13>=SI2EV 
HEHIX-J/1>-256 

fk;o~i. 0 

TRPE-l 

FIND F’fRSi LINE DESIRED 

RHI'JV0^--Rl=lWVl-I HflGVl*DELTflV 
Kl:=-1 ; OUTPUT LINE NO. COUNTER 

IF<RRWV0. QE. 0> GOTO 60 
DO 50 l^RRWVO. -1 
Kl^-Kl+1 

CRU.. 1PUTLNCK1> 

CONTI liUE 
RRI-iV0--0 

I P I'i l 1 7 VPE < TfiPE.. D I N < 0 > > 

DO 70 1=0.. RRWV0 
N08VTE=LINE<TRPE. DIN<0>, IFMT> 

IPXNDPVTE. EQ. 0) GOTO 930 
1F<RRW. EQ. 0> GOTO 225 

FIND SVNC MARKS 

ACCEPT "PRINT LINE HERDERS? <l-V/0-N>: ", IHD 

J =0 

1=29 

1 = 1+1 

I1=1PETCH<I.. DIN> 

I2=IFElCH<I+3, DIN> 

IFad. Gl. 200. AND, 12. LT. 50) GOTO 103 
IFCJ. GE. 5) GOTO 110 
rFCI.GH. N0BVTE-4> GOTO 130 
GOTO 80 

Jr-J+l 

SVNC<J)=I+3 
1 = 1+3 
GOTO 90 

TYPE "MULT SVNC MARKS FOUND AT : <SVNC< I >, 1=1, J) 
J--2 

CONTINUE 

IF<J, GT. 2) GOTO 110 
FIND MIDPOINT OF INPUT PICTURE 
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GO'IO 140, 150>.. J+1 

M]tV.SVNC<:i>+14S5 

IP<SVNC<1). GT. 1250> MID=SVNC<;i)-145S 
GO TO :160 

100 MID--<;SVNC';:i>-i-SVNC.<2>V'2. 0 

C INITIRLIiVP: CONSTRNTS 

160 S.--1 06666, ; SfiMFLIHG RRTE SRMPLES.''SEC0r4P 

Wr;,4ri0. ,• SCRM RRTE REVOLUTIONS/MIHUTE 

/CS. 0*P1*W> ,i SRMPLINQ COWSTRNT SRMPLES/RRDIRN 

COSQ:-0. 20791169 ,i ORBITRL INCLINRTION COSINE OF 7S DEGREES 
T-.vJ.44ti. J PERIOD OF EfiRTH ROTRTION MINUTES 
.! RRDIUS OF ERRTH kilometers 
DELPH l:---2, 0:t''::PI/W.'>*<l. 0/PERIOD+COSQ/T> ; RRDIRNS 
C0NSTd--<HEIGHT+R>7R CONSTRNT UNITLESS 

’ RLPHR~<RRWX1-MID>/K SCRN ANGLE OF REFERENCE POINT RRDIRNS 

Z-RSlNO;:ONSTli+iSINc:flLPHR>> ,• ZENITH ANGLE RRDIRNS 

PSI---Z-RLPHR ! GEOCENTRIC ANGLE OF REFERENCE POINT RRDIRNS 

THEIRi=RSIN<<SIN<.REFLRT*PI/180. >+COSO*SIN<PSI > >/COS<PSI > > i RRDIRNS 
DER.ETR:. 2. ©■^CPI/CW+I > )*SQRT<C0S<THETR)4!C0S<THETR)-C0S0^C0SQ> 
INTSL1-<1, 8ZDELPHI >-*.<2-RLPHfi-IMRGVl*DEBETfl*DELTRV) 

I N l X0r- 1 NTXI- I MRQX1*DELTRX 
LRSTX- < S I ZEX-l ^ '♦•DELTRX 
LINENO-'© 

LAST N0=^- < S I ZE V-1 ) DELTR V 
C WRITE OUT PARAMETERS 

WRlTECia.. 501> SIZEX, DELTRX.. SIZEV.. DELTflV.. HEIGHT, PERIOD, 

♦ RRWX:1., RRWVl, REFLR1.. IMRGXl, IMflQVl, INTX0, INTXi, 

=*< <THE‘IR>*=1S0/PI.’).. MID,. SVNC<1>, SVNC<2.'> 

001 FORMAT <1H1, 10X.. "SIZEX=", 17, 2X,. "DELTRX='‘, 17, 2X, ‘'SIZEV=", 

♦ 17, 2X, "DELTfiV=", 17, /Z, 10X, “HEIGHT^", FI©. 2, 2X, ■•PERI0D=“, 

Hi FJ.6. 2, 2X, ZZ, 10X, ••RRWX1=", 17, 2X, "RRWV1=", 17, 2X, "REFLflT='‘, FI©. 2, 

♦ ZZ, d OX, ‘MMfiGXl=“, 17, 2X, ''IMRQVl^", 17, 2X.. "INTX©=", F10. 2, 

♦ 2X, ”1^X1=", F10. 2, 2X, ••THETR=", F10. 2, ZZ, 10X, 

♦ "MID-", 17, 2X, "SVNC1=",I7, 2X, "SVNCZ^", 17, ZZZ> 

C PROCE.SS IMRGE 

170 K2==-l ; POINTS TO PIX IN OUTPUT LINE 

C0NST2=L INENO-t-DEBEI R 
DO 190 J=0, LRSTX,. DELTRX 
INTX---INTX0+J 

RNGLF> I NTX+DELPH I +C0NST2 

RRWX-- K4RTRN < S I N < ANGLE) Z< CONSTl-COS C ANGLE >>> +MI D 
C CHECK IF OUT OF BOUNDS IN X DIRECTION^ 

1F<RRWX. LT. 0. OR. RRWX. GT. NOBVTE-2) GOTO 180 
C USE LINEAR INTERPOLATION 
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IRW:=- 1FIX<RRWX> 
vi=a R m CH a RRWX. D IN > 
V'2!']Kb1CH<IRflNX+l. DIN> 
DRI.Ttt=-RflHX-IRRWX 
K2-K2+1 

I '.'HI, i-DELTR :> +V£-»=DELTR 

CRU S7 0RE<IVRL. K2. DRTOLIT) 

Gi:nc J90 ■ 

180 K2--K2+1 

CRI..1. S'l OREC©.. K£. Dfl10IJT> 

190 CONTINUE 

K 1 "-K 1+1 

CHI.L IPUTLN<Ki.^ 

IFCIHD. NE, 1> GOTO 198 
DO 192 1-0.. 27 

192 DR I OUT a ,'■ = I FETCH < I . D I N "> 

' WR1TF<12,. 601.T LINENO. <DRTOIJT<J>^ J=0. 

601 FORMHTCIH .. 15. 2'X,. £814) 

198 LINENO^-LINENO+DELTRV 

IKCK:1. GE, SIZEV-1) GOTO 990 

DO 200 1=1. DELTRV 

NOR't'TE- LINE CTREE.. DIN<0)> IFtIT) 

200 IFCNOB't'TE. EQ. 0) GOTO 210 

GDiO 170 
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C FILL H3TH 0 IF OUT OF BOUNDS IN V OR TAPE HRS FMT ERR 

210 IFCKl. GE. SISEV-1) GOTO 990 

DO £20 1=0,. SI2EX-1 
£20 CALL STORE<0. I, DRTOUT) 

DO £22 I=K1+1. SI2EV-1 
222: CRU.. IPUTLN<I.'> 

GOT 0 990 


I 


C 

225 

511 

£30 


RRW IHHGt 


£40 

£50 


£60 


I'lRlT t <10.. 511) 5IZEX.. DELTRX. SIZEV. DELTflV. RfiWXi. RRWVl 
FORI-iRTClH >616) 

KZ'-^ 

DO £00 ..I=RRWX1. RRWXl+CSI2EX-l>»=DELTflX> DELTRX 
K£. K 2+1 

IFC.J, LT. 0. OR. ..I. GT. NOBVTE-1) GOTO 240 
r..'Rl..= lFETCHC.J. DIN) 

CBI .1 ST ORE C I VRL. K2. DRTOUT ) 

GOT D £.50 

CRI..I ST ORE C 0.. K2.. DRTOUT) 

CONTINUE 

T't'PE Kl.. NOBVTE 

K 1 =K 1+1 

CRU,. IPUTLNCKl) 

IFCKl. GE. SIZEV-l) GOTO 990 
DO £.60 1=1. DELTflV 
NOt!VTE=LINE<TRPE. DINC0>. IFMT) 

IFCNOBVTE. EQ, 0) GOTO 210 
GOTO £30 
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('• ERROR MEIt'SHORS 

920 1VPF "TAPE FORMAT ERROR" 

QOfO 990 

940 TVPF "r-40 SVHC MARKS" 

GO'IO 991 

990 CALL WRBLK<2> 0.. HEAD, 1.. IER> 

CAI..L TIER<"WRHD'L IER> 

991 CUNI'JMUE 
CLOSE 1 
CL.OSE 2 
CLOSE 2 

CAM. RLSE<"MT0‘L IER> 

CALL TIER<"RLSEMT'MER> 
1000 STOP 

• ENI> 
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"MOAflSf-M-n" - SPECIFICflTlOM STRTEMENTS FClR PROGRAM NOflfl 

INI FfihR SIZEX. DELTRK.. SIZEV, DELTRV.. RAW.'-;!.. RflWVl, 

IMHGX1 . IMFiGVl. TAPE. SVNC:(,E«>.' BUF<236>,. RRWVO 
RFR!,. IH-1X0.. INTXl.. K 
IN1KGFR PRH 

IN’IF.CiFR DIN<0:?tIiO0.''j DflTOUT<Gi;SHi0> 

COHHONXCOLD/ DflTOU'1 

COl’Il'ini VPRRMSX SI2EH,. DEUTRH, SIZEV, DELTRV, RRWXl, RRWVl, 

» HE ] tiH r.. PERIOD, REFLRT, I MRQXl, I MRGV1> PI , RRW 

COHMON/NRMES/ UNIT<0;3:>, PRRT<0:S>, IPflRT<0:9>, NPRRT<0:6) 
INI fObk UNIT, PRRT, HERD 
COHHON/RRRflV/ HEflD<ei;255> 

DIMENSION R<120> 

• EQUIVALENCE HERD<16>> 
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SUBROLIIJNE IMPflR 


rPPT.S PRRRMETERS FOR NORfii 

IHCl.Ur'fc. "NORRSPEC" 

LOR/C.-RL NOTOK 

CRLI, a.H2LC0> 

TVPP "INPUT NORR PRRRMETERS<15>" 

RCi;:P;PT "SIZEK = ‘USISEK 
fiCCEPT "DELTRK = DELTRX 

ACCEPT "SI2EV = 'U SI2EV 
RCCTPI "DELTRV = DELTRV 

flCCEPl "RRNRl = RRNXl 
ACCEPT "RRWVl = 'C RRWVl 

ACCEPT "ENTER 8 FOR RAW IMAGE, 1 FOR CORRECTED 

IF CPRM. EQ. HC' GOTO 28 

IFCPHW. NE. 1> GOTO 28 

RCDEP'I "HEIGHT = ", HEIGHT 

ACCEPT "PERIOD = ", PERIOD 

RCCEP'I' "REFLRT = REFLRT 

RCCPP'I "IMRGXl = ", IMRQXl 

RCCEPT "IMRGVl = ", I MAG VI 

IF<NOTOK<0>> GOTO 10 

RE'ILIRN 

END 
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SUBkOUTlNE I OPEN 

C OBTRIN’^ IM'kf-CTORV AND FILE NAMES. 

CREATES AND OPENS FILE AND FILE. HD. 

0 ACCEPTS TAPE FILE NUMBER ArJD OPENS TAPE FILE 
C 

C TAPE OPEN ON CHAN 1 
O FILE. HD OPEN ON CHAN 2 
C FILE OPEN ON CHAN 3 

INCU,.iDF "NOAASPEC" 

LOOICAI. NOTOK 

CAI.I.. HITCH 

TVPE "ENTER DIRECTORV NAME" 

F'EAIXI:) . 521. ) NPARUG) 

FORl'iAT <S10> 

CALI. DIRCNPRRT. IER> 

CAI I,. PRaO< "PARTIN". I ER. 

C TEI..L. THE USER HHAT''S UP 

TVI-'E " " 

?20 PAU.SE I.OAD TAPE ON UNIT G. TVPE ANV KEV 
CAI .I.. JNIT C"MT0", 0. IER> 

GAIL FROCK" INIT'MER. *720) 

11 CAI..I.. Cl.H2L<0) 

C GFT THE TAPE FILE NO. 

LINITO:Oi--"MT" 

UNJT<l>--"0; *' 

30 TYPE "<15>".. "ENTER TAPE FILE NO. <0 IS 1ST? : " 

PF Al> <;t 1 . 1G00 > UN I T < 2> 

1000 F0FJ'1AT<S2> 

IF<NOTOK<0)>GOTO 11 
CALL. inOPtKl. UNIT. 1. IER> 

CALL ERGO<"TOPEN". lER. *30> 

C GET A DISK FILE NAME FROM USER 

730 DO 10 1-0.9 

IPART<1>=0 
10 CON'IINUE 

52 TYPE "<15>". "ENTER DISK FILE NAME. UP TO 10 CHARACTERS 
READCll. 1100 >PART COT 
1100 FORMAT <S10) 

rF<NOTOK<0T)GOTO 92 

C ADD THE , HD EXTENSION TO I PART 

CALL DOTHDCie. PARKO>. IPART<0>> 


71.0 • 

521 


CPffBTF OPEN THE . HL> FILE 

CRU.. CPlLl'KIPflRT. 2, IER;) 

IFcilER. EQ. l>QOTO 4 
IFOEP:. NE. 12>eOTD 5 
TVPh "HRME RLRERLVy IN USE" 

GOTO 92 

5 IFCIER. NE. 4>GOTO 6 

TVPE "TNVRLID NAtlE, ONLY fl-Z, 0-9 FIND ♦ PILLOWED" 
GOTO 92 

6 IFaPk. NE. 12.) GOTO 7 

TYPE "NO FILE WITH GIVEN NRME" 

GOT O 92 

7' CRI .!., EPJGO<"CRERTEHD". lER^ #730> 

4 CRI .L 0HEN<2. IPRRT.. 2, 512.. lER) 

Cfil . I.. E RQO < " HDOPEN " , I ER.. *720 > 

CRFHTE RND OPEN THE DRTR FILE 

CRI..I.. CFIlW<PRRT, 2.. lER,^ 

CRI..L FRGO<"CRERTEDRTR'S I ER, <:730> 

CRI 1. DPEN<2, PART, 2, SI2EX, lER) 

CALL TIER<'‘DTflOPN", IER> 

RETURN 

ENI.') 
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d.ci 


20 

2 t 0 


C 

40 

C 

50 

C 

60 

C 

70 

00 

9:10 

915 

920 

990 
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FUNCTION ITVPE< ICHRN, BUF> 

INIFOFK BUFCi)/ IR£C.f.‘l> 

COhlMON/NfiMES/ UNIKOiLO^ IDUM<26> 

INiFGFK UNIT 

1=0 

i=]+a 

L=rn4PFr'K53:40.. BUF,. ICHflN, IER> 

Cfll .1.. I) tP. < •' RDTPl ■’ , 1 ER ) 

IF<L, FQ. 45. OR. L. E0. 28-10. OR. L. EQ. 5240) GOTO 20 
IF*; I. FG!. 4) GOTO 91© 

GOTO 10 

IPFCa>=-L 
DO 2'0 1=2^4 

■ I RFC.': < I ) = I TRPEM < 52-10,. BUF. I CHflN.- I ER ) 

CRI.I.. 1 1ER<"RDTP2". lER) 

CONTINUE 

IPT=1 

IF<IRFC<1>. NE. 45. OR. IREC<2). NE. 45) GOTO 40 
IP! =2 
GOTO ^.0 

FORflRT 4 

IFTIRFCGIPT). NE. 5340) GOTO 50 

ITVPF=-1 

GOTO 990 

FORtIRT 1 

IFCIRFCCIPT). NE. 2840. OR. IREC<IPT+1). NE, 2840) GOTO 60 
I TVPF=d, 

GOTO 990 

FORflRT 2 . 

IF<rRFC<IPT). NE. 3840. OR. IREC< IPT-*-!). NE. 1500) GOTO 70 
ITVPE>-:< 

GOTO 990 

FORMRT 2 

IFCIRFCCIPT). NE. 45. OR. IREC<IPT+l). NE, 3840, OR. IREC<IPT+2>. 

♦ NE. 6) GOTO 80 - 

ITVPF=2 
GOTO 990 

IFGIPT. QE. 2) GOTO 920 

IPT=2 . < 

GOlO 40 

TVPF ''FRROR - AT LEAST 4 WRONG-LENGTH RECORDS” 

PFUJSE 
Cfll.l. BACK 

TVPF "ERROR - CANNOT HATCH TAPE TO NOa^ FORHAT " 

GOTO 915 

I Ct.1M---4O0Gi0K+5-IPT 

CRI I.. HTDI0<1, ICOft BUF, ISTS, IWCT> lER) 

CR(. I.. TJER<"TPBS",IER) 

RETURN 

END 
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INlHrGI-'R FUNCTION LINEaCHflN, BUF.. IFMT> 


READS (lliK l.INE FROM NOflfl TAPE 
LINE ^ HVIE COUNT IN LINfc READ 

INTFtiRR BUFcIl) 

LINFP^-U . 

L I NF1 - n RPEM 14009/ BUF < 1 > / I CHAN.. I ER ) 

CALI. FkGD<"RDTP".. lER/ T.iO> 

IF<UNE1. LE. e> GOTO 10 
GOTO <20,. 20/ 20,. £0>/ IFMT 

LINFP^ai APEM<14000-LINE1/ BUF< <LINEl+2>/2>, ICHAN/ IER> 
• CALL FRGD<''RDTP2"/ lER/ #10> 

IF<L TNK;2. LE. 6) GOTO 30 

LINE-L1NE1+LINE2 

IF<1. INF. LT. 2800 > GOTO 5 

RETURN 

END 


IV-24 


^3 




ORIGINAL PAGE IS 
OP POOR QUALITY 


INTKil-R KINCTIOr< I TfiPEM<NBVT. BUF, ICHRN. IER> 

INTt-GFR BUF<1> 

NWRI'3-G) 

Cfll L l'nDIO<ICHflN,‘ NHRD, BUF^ ISTS, IWCT. IER> 

CHI L r 1 ER < •' RDLN " , I ER > 

ITHPEH-2>^IWCT 

RETUR!^ 

ENI> 
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APPENDIX V, CASE STUDY 

IMPLEMENTATION OF THE REMOTE SENSING-AIDED PROCEDURE 
FOR ESTmTING SOLAR AND NET RADIATION AS APPLIED TO THE 
MIDDLE FORK OF THE FEATHER RIVER WATERSHED 

The stepwise procedure described in the procedural manual has 
been applied to a very diverse and large watershed in Plumas County 
in the Northern Sierra range of California known as the Middle Fork 
of the Feather-River Watershed (MFFR) . This watershed is approxi- 
mately 2900 km and is located at 39*^ 30' North to 40° 05' North 

latitude and 120° 05' West to 121° 25' West longitude . A Landsat 

color composite of the area with the watershed boundary on it is 
shown in Figure 1 . 

The areal distribution of all of the inputs to radiation models 
have been prepared and described models have been used to generate 
the final results. These results are net short-wave, net long-wave, 
and net radiation on the earth's surface for May 16, 1975. The 
input components consist of physiographic characteristics (elevation, 
slope, aspect, and albedo), together with climatic variables (cloud 
cover and temperature). Appendix II contains the results of 
this remote sensing-aided procedure for vegetation/terrain analysis 
as applied to the test watershed (MFFR). The vegetation map has 
been used to generate an albedo map of the watershed and is shown 
in Figure 2. 

Climatic data required for this remote sensing-aided system 
for estimating evapotranspiration was not fully available- within 
the MFFR watershed . The number of meteorological stations in this 
catchment is limited due to the sparsely populated nature of the 
area and the physical inaccessibility of a large portion of the 
watershed. Also, data pertaining to only a few types jaf climatic 
parameters have been taken at each station. There are only fif- 
teen ground meteorological stations within the Middle Fork of 
the Feather River watershed which collect climatic data. These 
stations along with their latitudes, longitudes and elevations are 
listed in Table 1. The data types collected at these stations con- 
sist mostly of precipitation and temperaturd information. Very 
few of the stations collect wind data or report on atmospheric 
condition with reference to degree of cloud cover. The main- climatic 
variable used in radiation models is the areal distribution of 
temperature. Two alternatives of interpolating temperature data 
based on ground meteorological stations (No. 1) and using NOAA-4 
satellite thermal data (No. 2) were Used to prepare the temperature 
data. Because of lack of adequate ground meteorological stations. 
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Figure V-1 
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Figure V-2. Areal distribution of albedo for the Middle Fork of the 

Feather River Watershed, in percent, based on vegetation/ 
terrain analysis 


Identification Key 


Color 

Albedo in 

White 

5 

Dark Blue 

9 

Purple Blue 

10 

Yellow 

11 

Light Green 

12 

Green 

13 

Dark Green 

14 

Light Blue 

15 

Orange 

18 

Red 

22 

Brown 

25 
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Table V-1. EXISTING GROUND METEOROLOGICAL STATIONS WITHIN AND SURROUNDING 

THE WATERSHED OF THE MIDDLE FORK OF THE FEATHER RIVER 


Station Name 


Brush Creek 

Bucks Creek 
Pump House 

Challenge Ranger 
Station 

N 

Forbes town 

Greenville Ranger 
Station 

Mohawk Ranger Station 

Oroville Ranger 
Station 

Plumas Creek State Park 

Porto la 

Quincy 

Sagehen Creek 

Sierra City 

Sierraville Ranger 
Station 

Vinton 

Woodleaf Oroleve 


Longitude,W LatitudeN 




39°41 


39°55 



X 


642.50 
643.02 

652.82 

649.53 

675.47 

704.17 

603.69 

697.23 

716.50 
675.94 

737.53 
702.74 

724.99 
740 . 64 
656.52 


Y 


4394.62 

4419.02 

4371.76 

4383.86 

4445.19 

4351.06 

4376.24 

4403.10 

4408.93 

4422.73 

4368.08 

4381.86 

4384.61 

4411.23 

4376.72 


ElevC-ft) 


3560 


Available Station 
Data for 1975* 
Preci " ' 
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alternative No. 1 is not recommended. However, a net radiation map 
based on the interpolated temperature values has been prepared and 
is shown at the end of this section. 

Using the procedure described in Step 8.b of the procedural 
manual, a temperature map of the study area based on the NOAA-4 
satellite thermal data and transferred to Landsat -based coordinates 
has been prepared. The date for this map was acquired as of 
9:30 AM PST^ May 16, 1975, and the temperature values range from 
+2°C to +4l”c. This temperature map is shown in Figure 3. 

Based on the procedure described in Appendix III, maps of 
elevation, slope, and aspect for the MFFR watershed have been 
prepared and are shown in Figures 4, 5, and 6, respectively. 

The elevation map is based on digitization of elevation contour 
isolines using USGS topographic series maps at a scale of 1:250,000. 

A map of net radiation for the study area (using the interpolated 
values of temperature from ground meteorological stations for May 16, 
1975) is shown in Figure 7. Another map of net radiation for 
the study area for May 16, 1975 (using temperature values based on 
NOAA-4 thermal data) is shown in Figure 8. 

Based on the evaluation of results it is recommended to use 
the NOAA-4 thermal data, unless adequate ground meteorological stations 
exist within the study area, for preparation of temperature input 
for radiation models;. Statistical summaries and percent distribu- 
tion of the data oVer the study area can also be calculated. As an 
example, percentage distribution of net radiation values in ly/day for 
May 16. 1975 over MFFR watershed (areal distribution is shown in 
Figure 7) is shown in Table 2. 
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Figure V-3. Temperature map of the Middle Fork of the Feather River 
Watershed for May 16, 1975, based on NOAA-4 satellite 
infrared data. 


Identification Key: 


Color 


TemperatuTe in 


o 


F 


Dark Blue 
Blue 
Purple 
Blue Green 
Green 

Yellow Green 
Yel low 
Yellow Red 
Red 

Lavender 


0-31 
32 - 45 
46 - 61 
62 - 68 
70 - 78 
79 - 84 
85 - 88 
90 - 92 
93 - 96 
97 - 108 
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Figure V-4. Areal distribution of the elevation zones for the Middle 
Fork of the Feather River Watershed, in feet (MSL) . 


Identification Key 


Color 

Elevation 

Red 

1200 

Orange 

1600 

Purple 

2000 

Brown 

2400 

Yellow 

2800 

Green 

3200 

Light Brown 

3600 

Light Yellow 

4000 

Yellow Green 

4400 

Dark Green 

4800 

Light Green 

5200 

Blue Green 

5600 

Dark Yellow 

6000 

Dark Blue 

6400 

Dark Purple 

6800 

Blue 

7200 

Purple 

7600 

Darker Purple 

8000 

Gray 

8200 


in Feet 
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Figure V-5. Slope map of the Middle Fork of the Feather River 
Watershed, in percent. 

Identification Key 


Color 

Slope in Pi 

Black 

Flat 

Red 

2-4 

Blue 

6-10 

Green 

12-20 

Yellow 

22-40 

Dark Green 

42-60 

Purple 

62-80 

Blue Green 

82-100 

White 

102-270 
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Figure V-6. Aspect map of the Middle Fork of the Feather River 
Watershed (N = North, S = South, E = East, W = West) 


Identification Key 


Color 


Direction 


Black 

Red 

Purple 
Orange 
Yellow 
Light Green 
Dark Green 
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Figure V-7. Net radiation map of the Middle Fork of the Feather River 
Watershed for May 16, 1975, using interpolated temperature 
values of ground meteorological stations as an input. 


Identification Key 


Color 

Net Radiati( 
in ly/ day 

Dark Green 

10 - 200 

Orange 

210 - 350 

Red 

360 - 400 

Blue 

410 - 430 

Yellow 

440 - 450 

Blue Green 

450 

Light Green 

470 - 500 

Purple 

510 - 560 
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Figure V- 


ORIGINAL Page IS 
OF POOR QUALITY 



Net radiation -nap of the Middle Fork of the Feather River 
Watershed for May 16, 1975, using temperature values based 
on NOAA-4 satellite infrared data as an input. 

Identification Key 


Color 


Net radiation in ly/day 


Dark Blue 
Blue-Green 
Yel low-Green 
Green 

Yellow-Brown 
Yellow 
Red- Brown 
Yellow-Red 
Lavender-Red 


0 - 250 
260 - 360 
370 - 410 
420 - 430 
440 
450 
460 
470 
480 


6i 
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Table V-2. 


it 


Percentage distribution of net solar radiation values 
in ly/day for May 16, 1975, for the Middle Fork of the 
Feather River Watershed 


Net Radiation 

10-200 

210-350 

360-400 

410-430 

440-450 

460 

470-500 

510-560 


Percent 

2.6 . 
6.2 
12.9 

19.5 

35.6 
12.2 
10.8 

0.2 
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V 

See Procedural Manual entitled, "Remote Sensing as an Aid in Watershed-Wide 
Estimation of Water Loss to the Atmosphere," pages 1-12 through 1-18 for 
references. 
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REMOTE SENSING AS AN AID IN WATERSHED -WIDE 
ESTIMATION OF WATER LOSS TO THE ATMOSPHERE 

INTRODUCTION 

Water vapor is the principal participant in the many energy 
exchanges taking place in the atmosphere. These energy exchanges 
are responsible for the weather phenomena which serve as important 
links connecting the various phases of the hydrological cycle. 

Quantification of the water flux mechanisms, viz. evaporation 
and evapotranspiration (which includes active plant water loss), 
is of importance in many scientific fields. They form one of the 
main components of the water budget, knowledge of which is indis- 
pensable for the solution of numerous water management problems. 

Reliable evaporation data are required for planning, designing 
and operating reservoirs, ponds, shipping canals, and irrigation 
and drainage systems. Evapotranspiration is especially important 
in arid zones where water must be used in the most efficient way 
possible. In addition, knowledge of the water requirement of 
crops depends partly on an accurate determination of the loss 
of water by evapotranspiration from cultivated fields. 

Evapotranspiration data are used for estimating irrigation 
requirements, rainfall disposition, safe yield of ground water 
basins, water yields from mountain watersheds, stream flow de- 
pletion in river basins, and in vegetation management for fiber 
production. 

In view of the overriding importance of evapotranspiration 
to accurate basin yield estimation and to other water consumpt- 
tive use determinations, any cost-effective, precise evapotrans- 
piration estimation system would be of significant utility to a | 

water manager. Remote sensing techniques may provide the key to I 

such a system's data requirements. Its timely, spatial, and 

relatively inexpensive nature give rise to accurate, location- F 

specific estimates of evapotranspiration, I 

I 

Evapotranspiration estimates are incorporated in current | 

water yield estimation procedures such as the joint Federal- i 

State River Forecasting Center model . It is therefore the purpose I 

of this section to lay the foundation for a remote sensing-aided | 

methodology for accurate and efficient watershed evapotranspiration j 

estimation. t 


Objective 


The objective of this procedural manual is to describe a 
stepwise remote sensing-aided procedure developed by our Remote 
Sensing Research Program (RSRP) personnel for the watershed-wide 
estimation of potential evapotranspiration from water supply 
basins of interest. The results of this study are designed to be 
readily implemented in hydrologic models for water yield forecasts. 

Definitions 

Evaporation may be defined as the transfer of water vapor from 
a non -vegetative surface on the earth into the atmosphere. Evapo- 
transpiration is the combined evaporation from all surfaces and the 
transpiration of plants. Except for the omission of a negligible 
amount of water used in the metabolic activities, evapotranspiration 
is the same as the ’’consumptive use” of the plants. The fact that 
the rate of evapotranspiration from a partially wet surface is 
greatly affected by the nature of the ground leads to the concept 
of potential evapotranspiration. Penman (1956) defines potential 
evapotranspiration as ’’the amount of water transpired in unit time 
by a short green corp, completely shading the ground, of uniform 
height and never short in water”. Pruitt (I960) designated the 
term ’’potential maximum evapotranspiration” to describe the situation 
when advected (horizontal movement of energy through the atmosphere 
surrounding the plant) energy is present. 

PROCEDURE 

This procedure is composed of a series of mathematical models 
for evapotranspiration (ET) estimation and techniques for providing 
the required input information to these models. T^e ET values as 
well as all the input values to the models are in areal distribu- 
tion (map) form. A complete description of a three level modeling 
technique and sampling design for providing the inputs to all three 
levels of models, as well as integrating the results, is given in 
Appendix I. The representative of level I models will be discussed 
in this manual. This representative model has been selected based 
on the comparison of the results of all level I models applied 
to a test watershed in Northern California. A modified version of 
Jensen and Haise model for daily prediction of potential evapotrans- 
piration has been selected for this procedure. Due to lack of 
funding, application of level II and level III evapotranspiration 
models has been postponed. 

1. Selection of study area and delineation of watershed 
boundary 

All the information for use in this procedure must be 
of a common register. In this case, the Landsat-base 



registration is recommended. Landsat multispectral scanner 
(MSS) data tapes of the watershed should be acquired and 
geometrically corrected. The watershed boundary should be 
determined and transferred onto the Landsat data base . The 
detailed description of the technique for geometric correc- 
tion of Landsat data and boundary transformation of the 
watershed is described on pages I-l through 1-8 (i.e., in 
Appendix I of the Procedural Manual for Net Radiation 
Estimation) . 

Selection of a model for evapotranspiration estimation 

The selected model is a modified version of the Jensen 
and Haise (1963) model. This model is defined as: 


ETp = (0.014T - 0.37) where 

ETp = daily potential evapotranspiration in inches 

T = mean daily temperature in °F 

R =' daily net radiation on the Earth's surface in 
inches of water equivalent 

Preparation of temperature map of the area as an input to 
the model 

Lack of adequate ground meteorological stations, in 
general, creates a problem in satisfying the temperature 
data requirements. We have chosen two alternatives for 
solving this problem: (1) computer assisted interpolation 
of existing temperature data collected in ground meteoro- 
logical stations; and (2) using the National Oceanic and 
Atmospheric Administration's NOAA-4 satellite Very High 
Resolution Radiometer (VHRR) data for generating the 
temperature map of the area. 

Collection and interpolation of temperature values 
based on ground meteorological stations data entail the 
following steps: 

a. Identification and selection of meteorological 
stations 

b. Registration of the data 

c. Recording and coding of the data 
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d. Temperature data interpolation 

e. Use of the data in simulation models 

f. Generation of final products 

The detailed description and the techniques for the 
above-mentioned steps are given on pages IV-1 through IV-4 
Ci.e. in Appendix IV of the Procedural Manual for Net 
Radiation Estimation) , 

Acquisition and processing of NOAA-4 satellite thermal 
data entail the following steps: 

a. Data acquisition 

b. Location of the study area on the raw thermal data 

c. Geometric correction of the data 

d. Data enhancement 

e. Transformation of voltage values to temperature 
values 

f. Transformation of temperature data to the desired 
coordinates system 

g. Preparation of final product 

The detailed description and the techniques for the 
above-mentioned steps and a brief description of the VHRR 
scanner are also given in Appendix IV (pages IV-4 through IV-26 
of this report) . 

4. Preparation of net radiation data as an input to the model 

The methodology for net radiation estimation is based 
on a series of mathematical models and a remote sensing-aided 
information system designed to provide the needed input data 
for the models. The detailed stepwise procedure for net 
radiation estimation is discussed in a separate manual titled 
"Remote Sensing as an Aid in Watershed-Wide Estimation of 
Solar and Net Radiation" appearing before the evapotranspira- 
tion manual. The results of net radiation models are site- 
and time-specific and consider the effects of climatic 
variability as well as terrain condition such as elevation, 
slope, and aspect. 
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5. Calculation, of potential evapotranspiration 

Upon the completion of temperature and net radiation 
maps, by using the equation described in Step 2 one may 
calculate the areal distribution of potential evapotrans- 
piration for the area of interest. T!ie computer program 
used for calculation of potential evapotranspiration (ETp) 
is known as MAPIT. MAPIT is a program developed by RSRP 
personnel to create, analyze, and display maps and map- 
related data. One may use his own technique for calcu- 
lation of evapotranspiration in map or any other form. 

An example of the results (potential evapotranspiration) 
as applied to the Middle Fork of the Feather Riv^r Water- 
shed, Northern California (approximately 2900 km'^) for 
May 16, 1975, is shown in Figure 1. (For a description 
of the study area, see Appendix V of the Manual on 
Radiation Estimation.) 

RSRP's film writing device known as the toage Ranged 
^tical Reproducer (IGOR) was used to generate the color- 
coded map shown in Figure 1 . The temperature map and net 
radiation map used as inputs for generation of this ETp data 
are shown in Figures 3 and 8, respectively, in Appendix V 
(case study) of the Manual for Radiation Estimation. 

NOAA-4 satellite thermal data has been used for preparation 
of a temperature map as well as being used as one of the 
inputs for preparation of the net radiation map. 
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Figure l Potential evapotranspiration (El'p) map of the Middle Fork 

of the Feather River Watershed for May 16, 1975 using 
temperature values based on NOAA-4 satellite data as an input, 
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THREE LEVEL EVAPOTRANSPIRATION ESTIMTION MODELS 
WITH REMOTE SENSING -AIDED SAMPLING PROCEDURE 
FOR PROVIDING REQUIRED INPUT DATA FOR MODELS 

The methodology for evapotranspiration (ET) estimation 
involves a multistage, multiphase sampling design employing 
three increasingly resolved levels of information to be used 
respectively in the three increasingly sophisticated levels 
of ET modeling. 

The sampling technique for providing all the required 
information for ET models (see the recommended models for 
their required inputs) is discussed at the beginning and is 
followed by the three levels of ET modeling. Data sources 
include the entire spectrum of environmental data gathering 
programs and information banks currently operating. Data 
types are those which are required by all three levels of 
ET models. 

1. Sampling technique 

a. Primary sampling units 

The first level is composed of Landsat and NOAA 
and topographic data. Vegetation-terrain and 
meteorological types of information are defined 
for a convenient base resolution element, in 
this case a small group of Landsat pixels. 

Based on this information associated with each 
base resolution element, an estimate of evapotrans- 
piration is made for that location. The appro- 
priate evapotranspiration equations, defined as 
level I models, perform adequately on this 
•'least resolved" information available from the 
first data level. 

After an estimate of evapotranspiration 
has been made for each base resolution element, 
the resolution elements (here 3x3, 4x4, or 5x5 
blocks of pixels) are grouped into primary sampling 
units . 


Secondary sampling units 

Within each primary sampling unit (PSU) 
selected, a series of secondary sampling units 
(SSU's) are defined and photographed at large 
scale with light aircraft. The photographic 
SSU data along with nearby snow course and 
ground station calibrated meteorological satellite 
data will provide the second level of information 
resolution. At the SSU stage of the sampling 
design, much more specific information concerning 
local vegetation canopy, soil, and local climatic 
conditions will be available as opposed to that 
obtainable from information at level I. Thus 
level II evapotranspiration models employing more 
data types and more refined data are used to 
generate evapotranspiration estimates for each 
PSU. 

Tertiary sampling units 

A sample of two of the SSU’s per PSU selected 
is then randomly (equal and/or unequal probability) 
chosen for further analysis on the ground. For 
each of these SSU’s selected, detailed ground 
measurements are made of vegetation canopy geometry, 
color, etc., as well as of soil and litter-organic 
debris conditions . The detailed data from this 
third level of information will then drive evapo- 
transpiration prediction models. Since estimates of 
evapotranspiration will be for the entire ground 
area of the SSU photo plot, this third stage 
unit (TSU) will actually comprise a double sample 
of the SSU’s. The sampling concepts and termi- 
nology in the preceding discussion are summarized 
in Figure I-l. 

Multistage equation 

A multistage estimation equation is being 
used to integrate the data from levels 1, II, and 
III for generating an estimate of watershed-wide 
evapotranspiration. This equation is shown in 
Figure 1-2 and the definitions of Its components 
are summarized in Table I-l. Based on this multi- 
stage equatioiV, watershed-wide estimation of 
ET is calculated by first using the ground-based 
evapotranspiration estimates to calibrate the 


Figure. I-l EVAPOTRANSPIRATION SAMPLE DESIGN TERMINOLOGY 

information 

Level Sample Unit Sample Design Level Comments 



STAGE I 


Multistage portion 
of sample design 


All PSU's in watershed 
have ^ evopotranspiration 
estimates made based 
on a summotion of 
evopotranspiration estimates 
for LANDS AT resolution 
elements occupying each 
PSU. A sample of the 
PSU's is selected for 
further sampling 



PHOTO SSU STAGE JL 



PHASE I 


Multiphase or diouble 
sample portion of 
sample design 


All SSU 's within a se- 
lected PSU have evapo- 
tronspiration estimates 
mode bosed on surface 
characteristics interpreted 
from the photography 



PHASE 2 Ground measurements 

resulting in ground- 
based evopotranspiration 
estimates are made for a 
sample of the photo SSU's 
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Multistage ET Estimation Equation 
and Corresponding Sample Frame 
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Table I-l Descript^^on of Terms Used in the Multistage ET 
Estimation Equation in figure 2. 

n = number of PSU's sampled 

= probability of selection of the ith PSU 

= number of SSU's selected within ith PSU 

JC {d| = the set of SSU's selected in each sampled PSU 

= area of the section within which SSU is located 

a = area of SSU 

Yjj = ground ET estimation for the jth SSU of the ith PSU 

p. , = probability of selection of jth SSU for ground 

measurement within ith PSU 

ET “ average evapotranspiration estimate for the entire 
watershed. 



SSU estimate derived from photo data. The cali- 
brated SSU estimates are then expanded to the PSU 
stage by utilizing the SSU selection weights 
developed earliei*. Finally, PSU evapotranspiration 
estimates can be expanded, each over the. appro- 
priate stratum, and then to the entire watershed 
by applying the PSU selection weights (proportion 
of evapotranspiration in the given PSU relative 
to all other PSU’s in the watershed) originally 
calculated. In this way, a cost-effective com- 
bination of an increasingly smaller sample of 
more precise and more extensive information levels 
can be utilized to give base-wide watershed esti- 
mates . 

e. Description of the procedure for using multistage 
equations 

First, the ground-based estimate of evapo- 
transpiration, for the SSU of the i^^ 

PSU is inflated by division with (the pro- 
portion of evapotranspiration estimated from 
level II data to occur in the SSU of the i^^ 
PSU) . After multiplication by the ratio of the 
average segment area (A;j^) to its corresponding 
plot area (a) , this inflated ground ET value 
represents an estimate of total evapotranspira- 
tion for the i'th psu. 

These PSU ET estimates are made for each 
ground measured SSU in the i^^ PSU. Ground 
visited SSU's in a given subsampled PSU, are 
defined to compose the set D. Summing over these 
ET estimates for each SSU in the set D and then 
dividing by n^^, the number of SSU’s in set D 
of the ith PSU, then gives an average estimate 
of PSU evapotranspiration. ITie steps just des- 
cribed are indicated within the bracketed ex- 
pression of the ET equation in Figure 1-2 . 

Once an average evapotranspiration estimate 
has been derived from level I and level II data, 
for each of the small number of PSU's subsampled, 
then each such value must be inflated to a water- 
shed-wide ET estimate. This is done, as shovm 
in Figure 1-2 by multiplying the bracketed quantity 
by There pi represents the proportion of 

evapotranspiration in the i^^ PSU based on level 
I information relative to the sum of evapotrans- 
piration from all PSU’s in the watershed. The 
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watershed -wide evapotranspiration estimates thus 
gained for each PSU selected for sub-sampling are 
then summed. This sum is next averaged by divid- 
ing by n, the total number of PSU's selected 
for subsamp.ling. The result is the total 
watershed-wide estimate of evapotranspiration. 

Sampling theory allows a derivation of an 
estimate of variance for the above watershed 
ET estimate. Multiplying the resulting standard 
error by an appropriate Student ’s-t value gives 
a confidence interval half -width (a statistical 
statement of estimate precision) to be associated 
with the estimate. 

The true value of evapotranspiration from the 
watershed is then expected to lie within plus 
or minus the confidence interval half -width of 
the ET estimate a given number of times out of 
a hundred (e.g. 90) with specified confidence 
probability (e.g. .95) . 

Recommended evapotranspiration models 

a. Recommended level I ET model 

A number of empirical and semi-empirical 
ET estimation models have been examined for 
this level. These models included Blaney- 
Criddle (1950), Hargreaves (1956), Turc (1961), 
and Jensen and Haise (1963) models. Based on 
the conq)arison of the results of the above- 
mentioned models with: 1) the pan evaporation 
data collected on the test site; and 2) expected 
values of ET for the test area, a modified ver- 
sion of the Jensen and Haise model was selected 
for the best performance. The test area was a 
catchment in the Sierra Nevada Mountains of northern ^ 
California, known as the Spanish Creek Watershed (475 km'^) . 

The modified Jensen and Haise model is defined 
as : 

ETp = (0.014 T - 0.37) R^ where: 

ETp = daily potential evapotranspiration in inches 
T = daily temperature in °F 

R^ = daily net radiation On the Earth's surface in Langleys 


1-7 




A complete stepwise procedure for this model is 
given in the main part of this manual. 

b. Recommended level II ET models 

The basis for the evapotranspiration models 
to be applied to the second level of information 
resolution is that of energy conservation, which 
is a proven law of thermodynamics. For this level, 
the energy-balance method will be combined with 
other methods for consideration of vegetation 
canopy effects and adveeted energy processes. 

The objective of the model applied at level II 
will be to capitalize on vegetation canopy, 
geometry- composition, and other surface data 
available from aerial photography to provide 
improved evapotranspiration estimates. The 
following models are being examined for application: 

Priestly and Taylor Model 

Priestly and Taylor used the equation: 

for estimation of evapotranspiration over the sur- 
faces where surface resistance was expected to be 
negligibly small. In this equation: 

S = slope of saturation curve versus temperature 
Y = psychrometric constant 
Rjj = net radiation 

G = heat flux into the soil 
a found to be 1 . 26 

This estimate of a is the overall mean of land 
and water surface . 

McNaughton and Black Model 

The best potential evapotranspiration relation- 
ship that was suggested by McNaughton and Black is 

Ph “ 0-17 GIL where: 

PE = potential evapotranspiration 
L = latent heat of water 

S = slope of saturation curve versus temperature 
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Y = psychrometric constant 

= net radiation 
N 

G = heat flux to the soil 
GIL = gross interception loss 

Modified Mcllroy and Slatyer 

The following equation was used by Mcllroy 
and Slatyer to estimate actual evapotranspi ration 
(E). 

E = § - (E . 9 -£ } — ^ where: 

s+Y h. 

1 

Y = psychrometric constant 

s = slope of saturation vapor pressure curve 
versus temperature 

h = transfer coefficient 

h. = molecular diffusion of water vapor through 
^ stomata 

= the canopy effect on the amount of actual 
^i evapotranspiration 

This equation will be modified to use some other 
terms for the canopy effect. 

Linacre Model 


A semi -empirical method of evapotranspiration 
was derived by Linacre (1967) through the combina- 
tion of Penman's method and solar radiation equations 
Linacre considered the forest canopy as in^jeding 
the loss of Water by means of two resistances: the 
canopy resistance, r , which is a function of canopy 
structures and their properties, and aerodynamic re- 
sistance, r , which is a function of the canopy- 
atmosphere interface. The combined equation is: 





A + Y + Y 


r 

c 

r 

a 


where: 
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difference between downward and upward 
radiation fluxes 

specific heat of air at constant pressure 


p = air density 
Y = psychrometric constant 
s = saturation deficit 

A = change of water vapor pressure with tem- 
perature in ambient conditions 


The variables to be estimated for input to the level 
II models are much more canopy and ground specific 
than those required for level I. In general, input 
variables to be applied in level II must be deter- 
mined more precisely and more frequently as well. 
Data requirements include temperature, humidity, 
and wind value profiles in the soil surf ace-vegeta- 
tion zone. Wind friction velocity, rainfall, ground 
cover, and topographic data may also be considered. 
These canopy data requirements will be fulfilled 
in large part through large scale photographic 
measurements collected form photo second stage 
sampling units and calibrated by corresponding 
ground data. Canopy- top meteorological data 
will be obtained as in level I, viz., through 
meteorological satellite data calibrated by 
ground station information and by microclimatic 
functions extrapolating satellite and ground data 
values to SSU's of interest. Vertical meteorolo- 
gical variable profiles within the canopy will 
be generated form canopy-top data according to 
canopy geometry functions. 


c. Recommended level III ET models 


The third information resolution level in the 
remote sensing-aided evapotranspiration estimation 
system will allow application of the practical as 
well as the most sophisticated models. The 
approach may be to select and develop those 
evapotranspiration estimation equations which are 
most rational and physical in terms of the actual 
processes involved. Based on the desired accuracy 
and need, a combination of empirical, energy balance 
and aerodynamic methods are recommended for exam- 
ination at this level. 
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For instance, this model (s) may be based 
on a combination of methods using radiation 
energy and temperature as primary variables 
(e.g. Jensen and Haise equation) and on 
methods considering vegetation canopy effects 
and aerodynamic resistance (e.g. Priestly and 
Taylor equation; McNaughton and Black equation) 
This new model (s) could be constructed so as 
to take best advantage of the data gained 
through remote sensing and ground sampling. 
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Figure 


Low altitude aerial photograph of the snowpack taken from 
a light aircraft using a 35 mm camera. The scale on the 
negative is approximately 1:30,000. The grid shown on the 
photograph represents an image sample unit (ISU) also 
located on the Landsat color composite imagery (See Figure 
1-4) . The image sample unit shown represents a snow cover 
class of 3 (See Table I-l) . . . . , 


Low altitude aerial photograph of the snowpack similar to 
that in Figure I-l. The image sample unit shown represents 
a snow cover class of 4 (See Table I-l). . 


August 13, 1972 Landsat color composite image showing the 
Feather River and the Spanish Creek Watersheds as they 
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April 4, 1973 gridded Landsat color composite image showing 
the Feather River and Spanish Creek Watersheds as they appear 
under maximum snowpack. The grids define the image sample 
units, each being approximately 400 hectares in size. . . . 

A portion of the April 4, 1973 Landsat color composite image 
showing the snowcover class interpretations (See Table I-l) 
on an image sample unit basis . . . . . , ... 
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This schematic diagram illustrates the photographic technique 

for constructing color composite images. A camera, color film, 

colored filters, LANDSAT imagery, a registration sheet and a 

light source are employed. Note that all three bands are 

placed separately on the registration sheet for copying. The 

band/filter combinations shown in the diagram above are used 

to produce simulated false color infrared composites. Different 

band/filter combinations will produce differently colored 

composites . . . . . . . . III-7 

Example exposure data for the construction of a simulated false- 

color infrared three band color composite using a conventional 

light table or a flashbox. LANDSAT MSS Band 4 is combined with 

a Wratten 47B filter. Band 5 with a Wratten 58 filter, and 

Band 7 with a Wratten 25 filter. Different exposure curves 

must be derived for different film speeds and/or different 

band/filter combinations as explained in the text.. . . . . . . . . III-8 
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REMOTE SENSING AS AN AID IN DETERMINING THE AREAL EXTENT OF SNOW 


INTRODUCTION 

The term "areal extent of snow" pertains to the surface area of the 
terrain in a specified watershed that is covered with snow at a given 
time or seasonal state. 

For many years the need for an, improved method of estimating snow 
areal extent has been recognized by hydrologists and others concerned 
with the prediction: of water runoff and the estimation of water yield, 
watershed-by-watershed. Among the methods showing greatest promise is 
the one dealt with in this Procedural Manual:, That method seeks to make 
maximum use of modern remote sensing techniques, including those which 
result in the acquisition of aerial and/or space photography of the 
watershed at suitably frequent intervals. Because of an advantage known 
as the "synoptic view," it is possible on space photography, and even 
on high-altitude aerial photography to cover vast watershed areas with 
only a very few frames of photography. As will presently be seen, this 
advantage can be of great importance to those wishing to use such 
photographs as an aid in determining the areal extent of snow. 

. I 

' ( ■ 

Of the several remote sensing-based techniques that have been 
developed for estimating snow areal extent, only two will be mentioned 
here since they appear to be both representative and highly promising. 
The first is based merely on the delineation of the apparent snowline 
on the photographs . The delineated portions of the photographs are then 
measured in such a way as to provide an estimate of the areal extent of 
the snow-covered area. This technique is described fully by Barnes and 
Bowley (1974). 

A second technique has been developed which utilizes a two-stage 
sampling scheme to arrive at an estimate of the snow-covered area. The 
procedure for employing this technique is described herein^ As will 
presently be seen, the technique employs as much ancillary (i.e., non- 
remote sensing-derived) information as possible to aid in the accuracy 
of the final estimate. Furthermore, the design of the model that is 
employed when this technique is used is such as to largely eliminate 
final estimate discrepancies , even when two or more photo interpreters 
are used as a team, since a statistical evaluation of each photo 
Interpreter's results is made. The method, overall, is sufficiently 
simple and inexpensive to permit it to be used in addition to, rather 
than instead of, more conventional methods, thereby providing additional 
confidence in the final estimate of snow areal extent on which 
operational decisions must be made, in any given instance, by the 
hydrologist or other resource manager. 
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When this second method is used by a competent image analyst, in 
addition to his determining the t otal area that is covered with snow, 
he also derives detailed locational information of great value, i.e., 
he determines the exact portions of the watershed that are snow-covered 
and differentiates them from the snow-free areas. This can be of great 
advantage in relation to his desire to make accurate estimations of 
water iyleld and water runoff rates because the rate of snow depletion 
for a given portion of the watershed obviously is a function of its slope, 
aspect, vegetative cover, and other locally-varlable factors, 

PROCEDURE 

The procedure described here for estimating areal extent of snow 
is based upon a two-stage analysis of remote sensing imagery. The 
first stage entails the use of Landsat* color composite imagery, while 
the second stage employs low-altitude, large-scale aerial photography. 

The following description provides a step-by-step guide to the use of 
this method in estimating the areal extent of snow: 

Step 1 . Prior to the start of the winter season (i.e., the season 
when the snowpack is to be measured) all of the projected dates when 
Landsat will be overflying the area of interest during the snow survey 
period should be systematically listed. In compiling this list one' 
should keep in mind that there is sidelap of about 20 to 30 percent 
(depending upon the latitude of the area) between the 115-mile swath 
that is covered by Landsat on one day and that which the same vehicle 
covers on the following day . Therefore, this factor needs to be 
considered as the list is being made of all possible dates for the 
acquisition of Landsat imagery of the area of interest. 

Step 2 . During the snow-survey period, on each date when a Landsat 
overflight of the area of interest is about to be made under what is 
predicted to be suitably cloud-free conditions, by prior arrangement 
steps should be taken to procure nearly simultaneous low-altitude 
aerial color photography of preselected flight lines, thereby acquiring 
the previously mentioned "second stage sample." An optimum scale for 
such photography is approximately 1/30,000, and 35mm negative color film 
is preferred. 

Step 3. The second stage photography procured in Step 2, above, 
should be promptly processed and from it standard "3R" prints (i.e., 
size 3-1/2" X 5") should be made. 


The term "Landsat" applies to various unmanned earth resources technology 
satellites (formerly called "ERTS") , each of which provides images of 
the surface of the Earth on an 18-day sun-synchronous cycle (0930 
local sun time, approximately) from an altitude of approximately 570 miles 
and covering a swath width of approximately 115 miles . 
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Step 4 . As promptly as possible after the corresponding Landsat 
overflight has been made, and again by prior arrangement (in this case 
normally with the EROS Data Center at Sioux Falls, South Dakota) black-- 
and-white MSS (multispectral scanner) transparencies should be obtained 
from bands 4, 5 ana 7 covering the entire watershed area of interest. 

The standard 9 ”x9" transparencies produced by the EROS Data Center are 
preferred . In addition, either at the time when the first of this snow- 
season Landsat photography is being acquired, or at some earlier date, 
an additional summertime (snow-free) set of Landsat imagery also should 
be acquired for comparative analysis, as indicated below (see Step 15) . 

Step 5 . In each instance, on the Landsat imagery that has been 
acquired the boundary of the watershed or other area of interest should 
be delineated. 

Step 6 . For each of the frames of Landsat imagery required in 
making the delineation of Step 5, a simulated color infrared color 
composite should be made, using the three bands (4, 5 and 7) to make a 
"triple exposure" on negative color film in the manner described in 
Appendix' II of this procedural manual. 

Step 7 . From the resulting negative color composites, reflection 
prints (rather than transparencies) should be made to a size of approximately 
8"xl0". In those instances where more than one Landsat frame is required 
to cover the area of intereist, the Landsat frames originally selected (and 
the corresponding areas copied and printed from them) should overlap each 
other sufficiently to provide complete coverage of the entire watershed 
or other area of interest. 

Step 8 . Either at the time when the first of the wintertime Landsat 
sets of imagery is obtained, or prior thereto, a set of acetate overlays 
each about 9Vx9" in size and gridded at various intervals into squares,^ 
should be prepared. In any given instance the size of the grid should 
be governed by time, accuracy, and cost constraints. Obviously, in 
most instances as grid size degreases, the accuracy, cost and time factors 
all increase. 


Step 9 . Through the use of geographic features that are readily 
located on Landsat imagery, each successive frame of such imagery should 
be Indexed so that the gridded acetate overlay can be registered in 
exactly the same way with respect to each of the frames that cover a 
given area . 

Step 10 . Similarly, the grids, as they appear on the Landsat frames 
of imagery, should be transferred to the low-altitude photos (i.e., the 
stage two photos taken at large scale with 35mm color film) and then 
labelled in such a manner that any grid unit (i.e., any image sample unit) 
can be easily located there and readily cross-referenced to the corresponding 
Landsat imagery. 


1. See: Table II-l Summary of Image Sample Unit (ISU) Optimization 
Data, page il-l 
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Step 11 . The indexed grids as they now appear on the low-altitude 
photography should be interpreted for areal extent of show. Generally 
four to six snowcover classes should be used, based on the proportion 
of the grid that appears to be snow-covered. An example using six 
classes could be: Class 1, 0-9 percent; Class 2, 10-25 percent; Class 3, 

25-50 percent; Class 4, 50-75 percent; Class 5, 75-90 percent; Class 6, 
90-100 percent. (It has been found that, due to the large amount of 
detail visible on the large scale prints, it is a fairly easy task for 
each grid unit to be placed in a specific class with a high degree of 
certainty.) 

Step 12 . This set of classified grid units should now be divided 
into a training and a testing set. The training set should be located 
on the gridded acetate overlays to the Landsat imagery and should include 
a wide range of snow appearance types covering the range of snowcover 
classes . 

Step 13 . The interpretation set-up should now be constructed. 

It should consist of a coincident image viewing device, such as a mirror- 
stereoscope or a transfer scope , the set of 8"xl0" Landsat summer color 
composite prints, for which the determination of snow areal extent is to 
be made, and the set of classified grid units, as seen on the low-altitude 
aerial photography, corresponding to the gridded winter color composite 
prints . 

Step 14 . As the detailed interpretation of the Landsat winter 
photos begins, one 8"xl0" summer Landsat print should be placed under 
the stereoscope (or transfer scope) along with the corresponding scene 
on the gridded winter Landsat print set. The prints should be adjusted 
so that conjugate images can be fused by the photo interpreter (i.e., 
so that all common features can be made to coincide with each other , one 
image viewed with each eye) , thus allowing the Interpreter to assess 
both the snow areal extent and the vegetation/terrain conditions in 
each grid unit. 

Step 15 . The interpreter should next engage in the task of training 
himself to recognize different snowcover condition classes within grid 
units based on: (a) underlying vegetation/terrain features (as seen on 

the summer Landsat print), (b) the snowcover visible on the gridded 
winter Landsat print, and (c) the appearance and subsequent classification 
of the grid! units as seen on the low-altitude aerial photography . By 
viewing through the coincident image device and observing how each grid 
unit in the training set appears during the winter and summer, and by 
then referring to the training set grid units on the low -altitude 
aerial photography, the interpreter can effectively complete this training 
task in very short order . 

Step 16 . Once trained, the interpreter should then proceed to 
classify each grid square (image sample unit) on the winter Landsat print 
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(In terms of the snow areal extent class to which It belongs) according 
to his synthesis of the Information available as he perceives it on the 
summer and winter Landsat prints. 

Step 17 . After all of the winter dates of imagery have been 
classified, the interpreter's initial interpretation results are ready 
for the testing phase. Therefore, at this point the results of the 
grid classification (as made on the low-altitude aerial photography 
testing set)i should be compared to the results obtained by the inter- 
preter in all applicable grid units. A ratio estimator statistic, as 
described in detail in Appendix I of this Procedural Manual should then 
be applied. This statistic will provide the adjustment calibration 
necessary to correct the Interpreter's initial interpretation results. 

Step 18 . Once the true value for the snow areal extent has been 
estimated using the population ratio estimator, two other statistical 
tests can be applied to clarify the final result. The confidence interval 
concept can be used to give an idea of the range of the final snow areal 
extent value. A Chi-square test can be applied to test the snowcover 
class width selection. Both of these tests are discussed in detail in 
Appendix I* 

Note: It is apparent from the foregoing that, even if this procedure 

is highly successful, its end product is merely an accurate determination 
of the areal extent of snow, together with an accurate delineation of 
each portion of the watershed area according to the snowcover class to 
which it belongs . This is an essential step leading to a prediction of 
water runoff and the estimation of water yield. The remaining steps are 
discussed in other procedural manuals. 



/' Appendix I 

^ J A SPECIFIC CASE STUDY 

The following case study describes uses that our RSRP personnel have 
made of the foregoing procedure for the Feather River Watershed, California. 

On April 4, 1973 the majority of the Feather River Watershed (located 
in the central Sierra Nevada Mountains, California) was imaged in a virtually 
cloud-free condition by the Landsat satellite. In order to document the 
existing ground conditions in greater detail than that provided by the 
satellite Imagery, a photographic mission was flown over the watershed on 
April 6, 1973., Four transects were flown using a motorized 35mm camera to 
acquire photography at a scale of approximately 1:30,000 on the negative 
(Figures ' 3 and 4) . ; 

Landsat imagery from August 13, 1972 was chosen for use as the 
vegetation/terrain base for the snowpack analysis procedures. The August 13 
and the April 4 Landsat images of MSS bands 4, 5, and 7 were photographically 
combined to produce simulated color Infrared enhancements (Figures^' 3 and 

4 ). . , 

The simulated color Infrared enhancement of the August 13, 1972 Landsat 
scene was photographically reproduced to give a 16” x 20” color print of 
approximately 1:250,000 scale. The April 4, 1973 Landsat color Infrared 
enhancements were enlarged to precisely the same scale as the August 13, 

1972 imagery and printed to give overlapping 8” x 10” prints. 

An acetate grid (each grid block equalling approximately 2,000 meters 
on a side at a scale of approximately 1:250,000) was attached to the 
vegetation/ terrain Landsat image. The grid blocks, termed image sample 
units (ISU), were transferred to the large scale photography where appli- 
cable. The image sample units on the large scale aerial photography were 
coded as shown in Table I-l. 


Table I-i 

Snowcover Condition Classes Used for Areal Extent of Snow Estimation 


Code Snowcover Class Midpoints 


1 No snow present within ISU 0 

2 0-20% of ISU covered by snow .10 

3 20-50% of ISU covered by snow .35 

4 50-98% of ISU covered by snow .74 

5 98-100% of ISU covered- by snow .99 


The interpreter then engaged in the training phase to familiarize him- 
self with the imagery and image classification technique. Once trained, the 
Interpreter proceeded to classify all the image sample units on the 8” x 10" 
Landsat color print for the winter dates. A more detailed description of this 
process will be found in the Procedural Manual. 
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FIGURE I-l 

Low altitude aerial photographs of the snowpack taken fron a light aircraft 
using a 35 mm camera. The scale of the negative is approximately 1’30 000 
The grid shown on the photograph represents an image sample unit fisui i 
located on the Landsat color composite imagery (see Figure 1-4).* The • ' ^ 
sample unit shown represents a snowcover class of 3 fSee Table I-l) 
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FIGURE IS 

August 13, 1972 Landsat color composite image showing the Feather River and 
the Spanish Creek Watersheds. 







FIGURE 1-4 

April 4, 1973 gridded Landsat color composite image showing the Feather River 
and Spanish Creek Watersheds. The grids define the image sample units, each 
being approximately 400 hectares in size. 
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The classified imagery was then subjected to statistical analysis. The 
following are the results from the April 4 Landsat areal extent of snow analysis. 



Snowcover (Condition) Class 

1 

2 

3 

4 

5 

Total number of ISU's 
classified per class 

403 

289 

214 

453 

859 


To derive the approximate number of hectares in each class, the image sample 
unit totals per class were multiplied by the number of hectares per image sample 
unit (400) with the results indicated below: 



Snowcover (Condition) Class 

1 

2 

3 

4 

5 

Hectares 

.161 X 
10^ 

: .116 X 
106 

.856 X 
105 

.181 X 
106 

.344 X 
10^ 


Then to establish the approximate amount of actual snow-covered area per 
class, the number of hectares/class were multiplied by the snowcover class 
midpoints (from Table I»-l) for each class, as Indicated below in the following 
tabular summary : 



Snowcover (Condition) Class 

1 

2 

3 

4 

5 

Snowcover condition 
class midpoint 

0 

.10 

.35 

.74, 

.99 .. 

Hectares of snow 
per class 

0 

.116x105 

.300x10^ 

.134x10^ 

.340x10^ 


By adding the "hectares of snow per class" values, the "total surface 
area of snow" in hectares was established. Thus for the April 4, 1973 Landsat 
daite the estimate of snow areal extent was found to be 515,772 hectares. 

This value, however, may be biased by any errors which occurred during the 
classification of the Landsat image sample units for snow areal extent. The 
uncorrected areal extent of snow estimate can be "corrected" by comparing the 
interpreter classification results with a series of preselected, preclassified 
image sample units, chosen for testing purposes from the large scale aerial 
photographs. 
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Interpretation results not falling in a given snowcover condition class 
for both the large scale aerial photographic and Landsat image data represent 
a misclasslflcatlon of image sample units by the interpreter (Table 1-2) . Thus 
it can be seen that 13 out of a possible 80 image sample units were mis- 
classified. This misclassif ication error can then be represented by a statis- 
tic, the population ratio estimator (Cochran, 1959), which will be used to correct 
the initial areal extent of snow estimate of 5.1582 x 10^ hectares derived from 
the Landsat image classification. 

The statistical approach is defined as follows; 

(Equation 1.0) Yj. = XR = the final, corrected areal extent of snow estimate 

N 

(Equation 1.1) where; X - ^ X. = the uncorrected Landsat areal extent of 

j=l ^ ' snow estimate 

given that; X. = the Landsat interpretation estimate of snow areal 
extent, per image sample unit, by snowcover class 

j = the index for all Landsat image sample units 

N = the total number of Landsat image sample units 
classified for snow areal extent 

(Equation 1.2) where; R = Y X = the population ratio estimator 






t 
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(Equation 1.3) given that; Y = 


(Equation 1,4) X = 


( y! yJ n = the average snow areal extent value 
1=1 / for the large scale photographic 

image sample unit estimates 

( yi n = average snow aceal extent value for 

i=l / the Landsat image sample unit 

estimates 



and that; n = the total number of image sample units used in the 
testing phase 

i - sampling index 

y^ “large scale photography snow areal extent estimate 
for image sample unit i 

= Landsat snow areal extent estimate for image sample 
unit i 
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Since the uncorrectki I-andsat areal extent of snow estimate! (X) has been 
calculated to be 5.158 x 4o^X hectares, the value for the population ratio 
estimator, R, must be de/iermined to solve for the final, corrected Landsat 
areal extent of i snow est^imate, Yj-i (from the equation Yj. = xS) . Table 1-2 
represents an expanded version of Table 2 designed to facilitate the calcu- 
lation of the population ratio estimator, R.- The table has been enlarged 
to include listings of snowcover condition class/image sample unit, specific 
large scale photography estimates of snow areal extent, and Landsat areal 
extent of snow estimates. 


= large scale aerial photography estimate of snow areal 
extent per image sample unit bj? snowcover condition 
class (in hectares) i 


snowcover condition class 

number of hectares 

midpoint J 

per image sample unit 


x^ = Landsat snow areal extent estimate per image sample 
unit by snowcover condition class (in hectares) 


= [snowcover condition, class j number of hectares ~\ 
I midpoint 11 image sample unltj 


fy = image sample unit interpretation frequencies 

u — index for table containing the results of the testing phase 
Additionally, 

n - the number of image sample units used in the testing phase, i.e. 


(Equation 1.5) 


- - Z fu = 8“ 

u=l 


where; u = index for Tables I “2 and 1-3. 

and k = the number of interpretation result blocks in Tables; 1-2 and 
1-3=11. 

S, the population ratio estimator, can now be calculated by solving for Y and 
X from the data presented in Table 1-3. Thus if 


(Equation 1.6) 
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(Equation 1.7) 

and 

r • 

n k 

Ti E ’ f u yu 

i=l U=1 

i 

OF POOR QUALITY 


then 

n 

Y- E yi + ^ 

i-l 

(Equation 1.3) 

(Equation 1.8) 


k 

= E ^u yu ♦ 

u=l 

i 

* 

(from Equation 1.6) 

. i. . ' ' 

(Equation 1.9) 


k k ^ 

= E (u yu - E'" 

U=1 U=1 

(from Equation 1.5) 


and 

_ n 

X = S Xi + n 
i=l 

(Equation 1.4) 

(Equation 2.0) 
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(from Equation 1.7) 
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Therefore (from Table 1-3 and equations 1.9 and 2.1), 


(0) (6)+(40) (l)+(40) (l0)+(296) (l)+(40) (2)+(140) (6)+(296) (2)-l-(140) (l)+(296) (12) 

+(296)(6)+(396)(35) 

6+1+10+1+2+6+2+1+12+6+53 

20784 


= 259.80 


and 




(0) (6)+(0) (l)+(40) (10)+(40) (1)+(140) (2)+(l40) (6)+(140) (2)+(296) (l)+(296) (12) 

: +396(6)+396f33) , 

6+1+10+1+2+6+2+1+12+6+33 

21132 

80 


264.15 



Consequently, the population ratio estimator (Equation 1,2) is calculated 


to be: 


R = Y + X 


= 259.80 7 264.15 


(Equation 1.2J 


= .9835 

I, i 

Finally, to solve for (the final, corrected areal extent of snow esti- 
mate), the population ratio estimate, R, is multiplied by the uncorrected Landsat 
areal extent of snow estimate, X (Equation 1.0). Thus, 


Yr = XR 


(Equation 1.2) 


" ( 515,772 ) (.9835) hectares 


“ 507,278 hectares 

Once the actual areal extent of snow has been estimated, Yr, it may be 
useful to know where the true value for the areal extent of snow lies. The con- 
fidence interval is useful for this purpose. The statistical approach is as 
follows, using the 95% level of confidence as an example: 

Yj. = the true value for the areal extent of snow 

The confidence interval around Y^. can be expressed as 
the probability that 




■r ; . 

(Equation 2.2) V-tn-1, 025^^^^^^^ " o25'^^^^ 

t L . 

(Equation 2.3) where V(Yr) [Y Yi + S " 2R £ x^ yi 


- sample variance 

n 2 ^ 9 

where Vy. = V f (Y ) 

1 ■f-', u 

1=1 1=1 i 

no k 9 

and V X = (X'^) 

1=1 1=1 

A summary of all the major results from the statistical analysis described 
previously for the April 4, 1973 snow analysis date is contained in Table VIII-4, 


I-ll 



IMiBl 








The final statistical test that was employed in this analysis was used 
to determine the suitability of the snowcover condition class widths as 
selected from step 11 in the procedure and Table 1-2. This test was to deter- 
mine if the values in the snowcover classes from the Landsat image data came 
from the same statistical probability distribution as the values in the snow- 
cover classes from the large scale aerial photographic data. If they came 
from the same distribution it could be expected that the areal extent of snow 
estimation procedure would provide good results. If these indicated that 
the two data sets came from different distributions, then new class widths 
should be selected and the areal extent of snow estimation procedure should 
be rerun. To perform such probability distribution "likeness tests," a 
Chi-square statistic was used: 

^^2 ^ (Oj-Ep ^ 

i=l ORIGINAL PAGE IS 

OF POOR QUALITY 

the observed values 
= the expected values 

k = the number of snowcover condition classes 
i = index. 

The values in the snowcover classes from the large scale photo data were 
designated as the expected values (El), since they were assumed to be "ground 
truth." The values in the snow cover classes from the Landsat image data were 
designated as the observed values (O-j ) . 

I 

I The results .of the Chi-square test for the April, 1973 statistical testing 
data (Table 1-5) are shown in Table 1-2, The results indicate' that the two 
data sets do come from the same statistical probability distributions; thus it can 
be assumed that the snowcover class widths used for the areal extent of snow 
analysis are adequate. 


(Equation 2.4) 
where 0^ = 
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Interpretation results of the testing phase, areal extent of snow 
estimation for the April 4, 1973 Landsat and April 6 ", 1973 large 
scale aerial photographic data. 


Large Scale AerlalJ Photographic Data 


Snowcover Condition Class 
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Key to Table X-2 


where: f^ = Image sample unit Interpretation 
frequencies 


u = Index for the table. 























Table 1-3 
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Interpretation results of the testing phase, areal extent ofi snow 
estimation for April 4, 1973 Landsat and April 6, 1973 large scale 
aerial photographic data. Also included is a listing of uncorrected 
large scale aerial photographic and Landsat image estimates of snow 
areal extent per image sample unit by snowcover class. 



where; Xu = Landsat snow areal extent estimate ? 

per image sample unit by snowcover [ 

condition class (in hectares) | 


Yu “ Large scale photography snow areal 

extent estimate per image sample ittvit 

I by snowcover condition class (iij hectares) 


f^ = Image sample unit interpretation frequency 
u =» Index for the table. 
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Table 1-4 

Summary of Statistical Results 
Areal Extent of Snow Estimation 
April 4, 1973 f 


Landsat estimate of the areal extent of snow (In hectares) 

X 

STS, 771 

Population ratio estimator 

R 

.9835 

Estimate of the true areal extent of snow (In hectares) 

/V 

Yr 

507.278 

Standard deviation of the areal extent of snow estimate 



(In hectare's) 

y-v(^r) 

12,659.29 

Confidence Interval (95%) around the true value for the 
areal extent of snow (In hectares) 


482 , 040sYr^ 532;5l7 

Total number of Image sample units Inventoried on the 
Landsat Imagery 

N 

2,218 

Average surface area of an Image sample unit (In hectares) 


400 

Total number of hectares inventoried 


887,200 

Total number of Image sample units sampled for the 
statistical testing procedure 

n 

80 
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Table 

Ghi-square Test for April 4i 1973 Data 

Null Hypothesis ; Ho; The observed values (0±) come from the same distribution 

as the expected values i (Hi) 

Alternative Hypothesis ; Hi; The observed values do not come from the same 

distribution as the expected values 

Significance Level : 5% 

Test Statistic under the Null Hypothesis; . 

1=1 ’ 

where k-1 = the degrees of freedom 
Hull Class Data Summary 


1:; ‘ ■ . •' 

Class 

Oi 

Ei 

(Oi-Ei)^ 

Ei 

r : 

1 

7 

6 

.1667 

'« ■*» i • 

2 

11 

13 

.3077 


3 

10 

7 

1.2857 


4 

13 

21 

3.0476 


5 

39 

33 

1.0909 


2 = 5.8992 
i-1 


Conclusion; The null hypothesis is accepted since the calculated 
value of 5.8992 is less than the table value of 9.49. 
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PROCEDURE FOR OPTIMAL IMAGE SAMPLE UNIT SIZE ESTIMATION 


II 


m if 


The operational application of the image sample unit areal extent of snow 
estimation technique is heavily dependent on the size of the image sample unit 
used. 

• It was expected that the smaller the ISU size the smaller the standard 
deviation (and variance) of the corrected areal extent of snow estimate. Con- 
versely, as the size of the ISU's decreased, it was suspected that inter- 
pretation time (also known as cost!) would ificrease. Iri an effort to give 
prospective users a starting point when deciding on whaf image sample unit 
size to use, a simple test was developed to find some optimal image sample 
unit size (i.e. minimized standard deviation and interpretation time). The 
test would compare interpretation time and standard deviation of the corrected 
areal extent of snow estimate for three different image sample unit sizes. 

One interpreter was used throughout the test to minimize - or at least stan- 
dardize - bias in the final results. 

Three image sample unit sizes (204, 459, and 815 'hectares each) 
were used on a single Landsat winter date. Interpreter training, classifi- 
cation, and testing were carried out as described in steps 15, 16 and 17. 

Tbe estimated values of the areal extent of snow along with their standard 
deviations and interpretation times for each image sample unit size, are summar- 
ized in Table' II-l, 

Table II -1 Summary of Image Sample Unit (ISU) Optimization Data. 
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Image Sample 
Unit Size (in 
hectares 

Standard Deviation 
of the estimate 
(in hectares) 

Interpretation 
Time (in hours) 

MV 

204 

9,981.36 

6.75 

> ?*•“ 

459 

16,864.54 

6.00 


815 

24,439.40 

5.50 

.f| 


As expected standard deviation decreases and interpretation time in- 
creases with image sample unit size. By plotting the data it was deter- 
mined that an image sample unit size of 365 hectares was the point 
where both standard deviation and interpretation were minimized as shown in 
Figure II-l 
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It should be noted that the image sample unit size of 365 hectares 
is to be considered a maximum size for 2 reasons: 

(1) The interpreter who analyzed the three image sample unit sizes 
was previously inexperienced in this type of interpretation. The 
largest ISU size was interpreted first and the smallest last. It can 
be expected that the interpretation time for each succeeding ISU size 
was reduced due to the experience gained on the operation of the in- 
terpretation system design. Thus, it seems reasonable that as more ex- 
perience is gained, the intepretation time for any ISU size would be 
reduced, lowering the'lnterpretation time versus ISU size" curve. Con- 
sequently, the value of 365 hectares for an "optional" ISO size 

could conceivably be reduced by increasing interpreter experience. 

(2) While the results illustrated in Figure II-l* .indicated an optimal 
ISU size for minimizing standard deviation and interpretation time, 

the image sample unit size may de decreased even further if more time 
can be spent in the interpretation phase. 
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Appendix III 


PROCEDURAL .MANUAL FOR THE PHOTOGRAPHIC CONSTRUCTION OF LANDSAT COLOR COMPOSITE IMAGERY 

INTRODUCTION 


The value of color composites has long been recognized as an aid to the 
interpretation of Lands at imagery. The technique described in this technical 
brief demonstrates how Landsat color composite imagery can be constructed 
from Landsat positive black-and-white transparencies conveniently and 
inexpensively,. 

COLOR COMPOSITE THEORY 

Typical color reversal and color negative films consist basically of 
three light sensitive layers. These layers react individually to blue, green, 
and red light, each layer corresponding to a particular dye, dependent upon 
film type, as indicated below: 


Light 

Color reversal film 

Color negative film 

Sensitivity 

Dye layers activated 

Dye layer activated 

Red 

Yellow and Magenta 

Cyan 

Green 

Yellow and Cyan 

Magenta 

Blue 

Magenta and Cyan 

Yellow 


If, for instance, color negative film is exposed to blue light, the 
blue sensitive layer of the film is activated while the other layers remain 
unaffected. This blue-sensitive layer corresponds to a final yellow dye 
image and hence the results will appear yellow on the negative. 

When color reversal film is exposed to blue light, the blue sensitive 
layer of the film is activated; however, during processing the two remaining 
sensitivity layers corresponding to the magenta dye (green sensitive) and 
the cyan dye (red sensitive) are either chemically or mechanically exposed. 
The originally exposed blue sensitive layer becomes non-functional and the 
resulting image appears blue due to the combination of magenta and cyan 
dyes. The reactions of color reversal and negative film to red and green 
light can be explained similarly since red is produced by a combination of 
yellow and magenta dyes and green is produced by a combination of cyan and 
yellow dyes. 





The photographic construction of Landsat color composites relies on the 
three-layer sensitivity of color film and on the different densities of any 
given feature as it, appears on the various Landsat black-and-white bands. 

Esch of the different Landsat black-and-white bands chosen to make the color 
composite is photographically copied onto the same sheet (or frame) of 
color film, each baiid through a different colored filter (in essence a 
multiple exposure). ' Figure III-l illustrates the set-up used to produce 
Landsat color composites. Depending on which filter is used with each band of 
imagery different color renditions miay be produced, hopefully enhancing 
; different features. 

MATERIALS AND EQUIPMENT 
j LANDSAT BLACK-AND-WHITE IMAGERY 

‘ The construction of Landsat color composites commonly uses Multispectral 

! Scanner (MSS) imagery of bands 4, 5, and 7, usually in the 9h inch 

positive transparency form (although the 70 millimeter positive 
I transparencies may also be used). MSS Band 6 is also used for some 

I composites as well as available Return Beam Videcon (RBV) positive 

imagery. Furthermore, just as a color composite can be formed by 
combiniiig bands of imagery acquired on a single date, so it also can 
be formed by combining several dates of imagery acquired in a single 

band, or any combination of multiband/multidate images. 

! . , ... 

FILTERS 

■ Eastmak Kodak Wratten filters are suggested for use in the construction 

I of color composites due to their optical quality and the large variety 

of colors and sizes available. To produce a simulate^ false-color in- 
frared color composite image (Figures 1-3 and 1-4) , Wratten fil- 
ter number 25 (tricolor red), 47B (tricolor blue), and 58 (tricolor green) 
can be used. 

i 

FILTER RELATED ACCESSORIES 

1 Eastmak Kodak filter frames provide a convenient way to mount the 

i Wratten filters. The filters so mounted can be easily inserted into 

' a commercially available filter holder which mounts, via adapters, to 

the filter- threaded receptacle in the front of most enlarging lenses. 

{ ■ ■ ... 

PHOTOGRAPHIC LIGHT SOURCE 

! The light source is used to provide the light which will be transmitted 

I through the Landsat bands during the exposure process. Basically, the 

light source may be of two common types; incandescent or flourescent 
and electronic strobe. Incandescent or flourescent illumination in 
the form of a "light table" may be the most convenient form to use 
i initially. One of the problems encountered with using a light table 

! is its relatively low light output, necessitating long exposures times. 

One way of countering this is by using an electronic strobe light box 
i which could give decidely greater light output and thus shorter exposures. 
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REGISTRATION SHEET - 

The registration sheet is used to register the individual Landsat bands 
with respect to one another. It also serves as a convenient base for 
introducing grids, watershed boundaries or other data directly on the 
color composite. The material used for the registration base should be 
translucent (frosted preferably), dimensionally stable and have a surface 
amenable to scribing. 

SCRIBING TOOLS 

The scribing tools consist of a scriber and a metal straight edge. The 
scriber can be constructed by inserting a drafting compass needle in an 
engineering drafting pencil while the metal straight edge may siraly be 
a metal engineering scale. The scribing tools are used to transfer the 
Landsat registration marks onto the registration sheet. 

CAMERA SYSTEM 

A large variety of cameras and related accessories ranging in format 
from 35 millimeter to much larger image sizes can be used to construct 
Landsat color composites. The only requirement is that the camera be 
capable of making accurate multiple exposures on a single sheet or 
frame of film. The most convenient camera used to develop this system 
was found to be a type of copy camera which utilizes a leaf shutter 
system allowing as many exposures as desired on a single sheet or 
frame of film. A camera system of this type is desirable due to its 
ability to accept film formats from 35 millimeters to 4x5 inches. 

The camera lens used for making color composites should be of very high 
quality. Since the construction of color composites utilizes most of 
the light-sensitive spectral range of conventional color film, the lenses 
used should be color corrected. The lenses should be of very high 
uniform edge-to-edge resolution and contrast to assure that no 
information is lost throughout the image area. 

EXPOSURE METER 

An exposure meter should be used to determine hov? much exposure should 
be given to each of the Landsat black-and-white images in making the 
color composite. The meter used should have a locking needle, a 
wide sensitivity range and a narrow acceptance angle so that selective 
area measurements can be made. 

FILM TYPES 

If reflection prints are to be made of the color composite, it is most 
convenient to use a fine grain negative color film such as Kodak 
Vericolor II Type S. If positive color transparencies are desired (e.g. 
for projection purposes) the use of either Kodachrome or Ektachrorae film 
is appropriate, depending on the film format used. 
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PROCEDURE 


The following step-by-step description will detail the procedure used 

to construct a Landsat color composite. 

Step 1 . Tape one band of a Landsat black-and-white set of imagery for 

a specific date to a light table in proper viewing position. 

Step 2. Tape a sheet of the registration base material over the Landsat 

image so that the four registration marks C+) iti the corners of the 
Landsat image are completely covered by the overlying sheet. 

, 

Step 3. Transfer the registration marks to the registration base material 

using the scriber and the metal straight edge. The scribed line should 
be narrow in width and extend beyond the actual Landsat registration 
marks to facilitate the actual registration process. 

Step 4. The registration sheet should be labeled with the latitude /longitude 

and date of the Landsat imagery from which it was made . This information 
is placed on the same edge as observed on the Landsat image below so that 
the Landsat black-and-white bands may be oriented correctly during the 
registration process. 

Step 5. The completed registration sheet is removed and placed in proper 

position on the light table orflash box and taped down. 

Step 6. One of the Landsat bands to be used for making the composite is 

taped down in register over the registration sheet. The registration 
marks and the information line are used to register and orient the two, 
relative to one another. 

Step 7. The camera system is focused upon this registered image and 

locked into position. 

Step 8 . An exposure measurement of the Landsat hand is then made on the 

ground glass focusing screen of the camera with the lens set at a fixed 
aperture, usually maximum. This reading is recorded. The other Landsat 
bands to be used are then placed one at.,ia time on the registration sheet 
with the image projected onto the focusing screen and the exposure 
measurement for each is recorded again at the same fixed aperture. 

Step 9. The exposure measurements for the Landsat bands are averaged 

and refered to exposure graphs (which are explained in the next section). 

Step 10. The camera/Iens exposure controls are then set according to the 

exposure graphs . 

Step 11. The film on which the composite is to be made is then loaded in 
the camera. 
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12. One at a time, the Lands *1; bands are placed on the registration 

sheet and exposed through the appropriate filter. The film is not 
advanced or removed until all exposures have been made. 


composite, 


The film is then processed normally yielding the final color 
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EXPOSURE DETERMINATION 

'In determining t^e proper exposure; to give each of the Landsat black-and- 
white images the separate bands are assumed to be matched separation positives. 
The actual deviation from this assumption is generally so slight that it 
may be ignored in most if not all cases. 

Since no two light source systems developed to construct Landsat color 
composite imagery will have equal light outputs, the graph used to determine 
proper exposure must be empirically established. The following describes 
the determination of the exposure graph necessary for proper exposure calcula- 
tion. Since each exposure graph can be calculated for one particular color 
composite type Ci*e., band/filter combinatioh) only, the following description 
will deal with the construction of a simulated false color infrared composite. 

To make this particular rendition Landsat band 7 is exposed through a 
Wratten 25 (Red) filter, band 5 through a Wratten 58 (Green) filter and band 
4 through a Wratten 47B (Blue) filter. 

To empirically determine the exposure graphs for the various band/filter 
combinations, a series of color composites are constructed at different 
overall exposures. Each series of color composites has exposure values measured 
for each band as per step 8 under "procedures'*. The values measured are 
averaged for each set and then plotted versus the actual overall exposure 
found to give an acceptable color composite. All light measurement values are 
taken with the lens set at a fixed aperture. This aperture will always be used 
for the through-the-lens exposure measurements of each band. This is 
necessary to provide a consistent index system for reference. to the exposure 
graphs . 

The averaged exposure measurements taken from the individual Landsat 
bands are plotted on the x-axis versus Exposure Values (EV) on the y-axis . 

EV refers to equal light intensities; hence equal combinations of shutter 
speed and aperture which yield the same exposure to the film. For example, 

EV 13 equals 1/30 second at f/15, .1/60 second at f/11, 1/125 second at f/8, 

1/250 second at f/5.6, etc. Figure II 1-2 illustrates the appearance of such 
an exposure graph. 

Due to the different optical densities of the filters used to construct 
a color composite, different exposures are needed for each band/filter 
combination. These different densities can be refered to as filter factors , 

The filter factor for a given filter changes with the film type as well as 
the light source; consequently filter factors must be calculated for each 
film type/light source combination. One practical way of calculating specific 
filter factors is to photograph a transparent gray scale, through eacl’ of the 
filters for which the filter factor is to be determined. These will produce 
individual exposures through the filters onto color film, either positive or 
negative. Several exposures are taken through each filter, advancing or 
changing the film after each exposure. All exposures must be recorded for 
the final determination of the filter factors. For each filter an exposure 
is found which best reproduces the gray scale (i.e. all gray levels are 
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clearly distinguishable from each other) by examining the color film. In all 
cases each individual gray level should have the same optical density for all 
of the filter types used. This is best done through the use of a densitometer 
but can be done visually if necessary. A comparison of the exposures used for 
the various filters to obtain equal densities gray level-by-gray level, will 
give the filter factors for each filter, These filter factors may then be 
translated into EV’s to give the relative exposure differences for the filters 
used. 

To establish the exposure graphs for the various band/filter combinations 
a series of color composites are constructed, all at different overall 
exposures . Each composite so constructed must have the filter factors (in 
terms of individual exposure differences) included. If, for example, for 
a specific film/light source combination, it is determined that a Wratten 
47B filter requires I EV more exposure than a Wratten 58 filter, (which 
requires itself H EV more exposure than a Wratten 25 filter to obtain equal 
density levels) then a series of exposures can be made using these differences 
For the first composite the overall exposure might be EV13 through the Wratten 
47B filter, EV 14 through the Wratten 58 filter, and EV 14.5 through the 
Wratten 25 filter. For the second composite overall exposure might be EV14 
through the Wratten 47B filter, EV 15 through the Wratten 58 filter, and EV 
15.5 through the Wratten 25 filter and so on. Since the individual filters 
require different but equal relative exposure differences maintained during 
the constiuction of the composites, color balance is kept stable, but overall 
composite density is not. Thus the best composite density from this series 
of different overall exposures can be determined. By plotting the exposures 
used for each filter in the best composite constructed, one group of points 
may be plotted on the exposure graph for this particular band/filter 
combination. Once two groups of points, each representing a good color 
composite overall exposure for a specific band/filter combination have been 
plotted, the graph can be completed. 

CONCLUSIONS 

Several critical points are worth mentioning again to assure that the 
color composites are correctly constructed: 

1. During the actual photographic construction, several exposures 
are made on the same sheet or frame of film; the film is not 
removed or advanced during these exposures. 

2. All exposure measurement values of the various bands are made 
through the camera/lens assembly on the ground glass using a 
fixed aperture. This aperture may change depending upon the 
actual responses of the various bands during the copying process . 

3. The individual exposure measurements of the bands to be used are 
mathematically averaged and then referred to the exposure graph 
for the exposure setting of each band/filtei combination. 
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CAMERA SYSTEM 



Figure III-l This schematic diagram illustrates the photographic technique for 
constructing color composite images. A camera, color film, colored filters, 
LANDSAT imagery, a registration sheet and a light source are employed. Note 
that all three bands are placed separately on the registration sheet for 
copying. The band/filter combinations shown in the diagram above are used to 
produce simulated false color infrared composites. Different band/filter 
combinations will produce differently colored composites. 
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Figure III-2 Example exposure data for the construction of a simulated false-color 
infrared three band color composite using a conventional light table or a 
flashbox, LANDSAT MSS Band 4 is combined with a Wratten 47B filter. Band 5 
with a Wratten 58 filter, and Band 7 with a Wratten 25 filter. Different 
exposure curves must be derived for different film speeds and/or different 
band/filter combinations, as explained in the text. 


*Exposure values CEV) are those unique combinations of shutter speed and 
aperture which yield the same exposure. For instance EV 13 represents the 
combinations 1/125 second at f/8, 1/60 second at f/11, 1/30 second at f/16, etc. 
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Procedural Manual 


REMOTE SENSING AS AN AID IN DETERMINING 
THE WATER CONTENT OF SNOW 

1,0 Introduction 

This procedural manual describes the use of Landsat data in 
combiTiation with conventional ground snow course data to provide 
an estimate of watershed or subbasin snow water content. The tech- 
nique is designed to be implementable in operational water runoff 
forecasting models incorporating, as one variable, an estimate of 
snow water content (expressed as a depth of water or as a dimen- 
sionless index) . Developed and tested in the Northern Sierra 
Nevada Range in California (Thomas and Sharp 1975), the method may 
find its greatest utility in mountainous situations characterized 
by significant winter snowpack buildup. 0nly a small initial 
capital investment is necessary and personnel requirements need 
not be substantial. Tlie technique is designed to complement 
current snow measurement methods by providing spatial information on 
snow water content. Consequently it permits more accurate estimates 
to be made on a basin or subbasin basis. However, successful 
application of the method does require relatively cloud free weather 
following winter storms during which sufficient Landsat coverage 
can be obtained of the watershed (s) of interest. 

Normally, snow water content estimates are obtained directly 
from ground-based snow course or snow sensor measurements. The 
procedure described herein introduces a stratified double sampling 
approach that relates the ground-based estimates to snow areal 
extent data gathered from Landsat imagery. The resulting relations- 
sh ip s enable low-cost remotely sensed data to characterize statis- 
tically the spatial and temporal variability of specific snow 
depletion environments sampled with ground snow course data. In 
this manner satellite data can be used to determine the weight 
assigned to a particular snow course measurement and also to provide 
more frequent assessment of snow water content. 

For determining snow areal extent, this technique utilizes 
the Landsat-based procedure described in our Snow Areal Extent 
Procedural Manual as the remote sensing input. However, non- 
Landsat remotely sensed data types also could provide useful 
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information to characterize the spatial variability, watershed-wide 
with respect to snow water content. For example, meteorological 
satellite data could be used with the Sncw Areal Extent Procedure 
to provide more frequent (daily) information, although of lower 
spatial resolution. The technique can also be refined by the user, 
if desired, to include machine processing of the satellite data. 

2.0 General Approach 

The following provides an overview of the remote sensing-aided 
snow water content estimation procedure. 

Sample Design and Measurement 

A stratified double sample method is used to develop a basin- 
wide estimate of snow water content. Under this approach, snow 
water content information for the whole watershed, as obtained in- 
expensively on a sample unit basis from Landsat data, is combined 
with that gained from a much smaller and more expensive sample of 
ground-based measurements at snow courses (see Figure 1) . The 
result is a basin-wide estimate of snow water content based on Land- 
sat data calibrated by regression on ground snow course date , Since 
much of the watershed snow water content variation is accounted for 
by information gained from the Landsat sample stage, an overall 
estimate of basin snow water content is possible at more precise 
levels (i. e., lower sampling error) than available for the same 
cost from conventional snow course data alone. 

The sequential sampling and measurement process proceeds by 
first locating a sample grid over the watershed. Snow areal extent 
estimates are quickly made for each sample unit by manual techniques 
(see Snow Areal Extent Procedural Manual) for the previous snowpack 
buildup dates and then for the specific forecast date. The snow 
areal extent data is then combined by a linear equation to generate 
an index parameter that is correlated with snow water content in- 
formation specific to the forecast date for each sample unit. 

This linear model is designed to reflect the relationship between 
snow areal extent and snow depletion behavior, and is specific to 
the Watershed being studied. Some users may choose to develop more 
complex, physically realistic areal extent-to-water content trans- 
formations . 

By specifying the precision (sampling error) and level of con- 
fidence desired in the basin-wide snow water content estimate and 
by considering measurement costs in relation to the available budget, 
one is able to calculate the necessary ground subsample size. Ground 
snow water content measurements are then allocated to Landsat-based 
snow water content-index classes (strata) according to weighted ran- 
dom stratified sampling procedures. 


STRATIFIED MULTIDATE LANDSAT DATA PLANE 
CALIBRATED BY SNOW COURSE MEASUREMENTS 
FOR WATERSHED SNOW WATER CONTENT ESTIMATION 
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Regression relationships are developed between the Landsat 
snow water content index data and the ground snow water content 
measurements. These equations are then used to correct all Land- 
sat-based data by ground values of snow water content. The ground 
corrected values of Landsat-based snow water content information in 
each stratum or class are added to give a total basin-wide estimate 
of snow water content, together with an associated confidence interval 
around the estimate. 

Results and Their Applications 

The Landsat-aided snow water content estimation procedure is 
designed to geherate an estimate of total watershed snow water con- 
tent and an associated statement of precision for a given forecast 
period. The estimate may then be related by regression equations 
directly to basin water yield for a given period. Alternatively, 
the statement may be used as another predictor variable in current 
snow survey or river forecast equations. Grid overlays, when placed 
either on watershed maps or directly on Landsat imagery, give a 
location- specific estimated snow water content index for each sample 
unit. Such information can be used to produce improved hydrologic 
modeling procedures incorporating this spatial data. 

Performance Criteria 


Performance criteria for this procedure consist of 1) the 
precision (sampling error) and accuracy of the estimate of snow water 
content over the watershed versus the cost of making the estimate; 

2) the timeliness of the estimate; and 3) the value of any in-place 
map information regarding snow water content that can be produced 
through use of the procedure. 

3.0 Detailed Description of the Method 

The stepwise procedure for estimating snow water content with 
the aid of remote sensing is described below: 

Step 1 : An overall plan is first developed which will facilitate 

adoption of the remote sensing-aided snow water content 
estimation procedure by the user. This implementation 
plan should consider 1) available budget; 2) requirements 
for either training or obtaining image interpreter (s) ; 

3) type of products desired, e.g. , watershed and/or 
sub-watershed estimates or in-place snow water content 
maps; 4) performance requirements in terms of estimate 
precision and satellite image acquisition- to-forecast 
turnaround time requirements; 5) interface of the snow 
water content procedure with current user forecasting 
methods; and 6) startup equipment (stereoscopy, photo- 
graphic laboratory facility) and labor requirements. 
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Step 2 : Color composites are prepared from Landsat imagery 

and/or from other satellite imagery types (e.g., 
weather service) , and the appropriate image sample 
unit (ISU) grid placed over each watershed of interest. 
In the performance of this step, black-and-white 
Landsat transparencies are obtained (e.g., from the 
EROS Data Center) and transformed into a simulated 
infrared color composite by a sample photographic 
procedure described in the Snow Areal Extent Pro- 
cedural Manual. In the color- combining process, an 
ISU grid is randomly placed over a base-date image to 
cover the watershed of interest. The base date Land- 
sat transparencies should be a summer image clearly 
depicting watershed geomorpho logical and vegetative 
features obtained before the first winter in which 
this technique is applied. Grids on all succeeding 
Landsat transparencies subjected to the color-com- 
bination process should be aligned with that on the 
base date image. In addition, watershed boundaries 
as plotted on 1: 250,000-scale US Geological Survey 
topographic map sheets should be transferred to the 
color composite images via optical rectification de- 
vices (if available) or by ocular inspection. 

A complete description of the boundary and sample 
grid overlay is given in the snow Areal Extent Pro- 
cedural Manual. Square image sample units each 400 
hectares in size, are recommended. This size is large 
enough to stabilize variance and small enough to give 
location-specific snow water content-related informa- 
tion at reasonable cost. The user may choose to 
modify these image sample unit dimensions as experience 
is gained with the estimation procedure. 

Step 5 : The snow areal extent is then estimated for each 

Landsat ISU for previous season (s) or current season 
snow buildup dates. Each ISU should be interpreted 
manually to determine its average snow areal extent 
cover class according to the snow environment-specific 
technique described in the Snow Areal Extent Procedural 
Manual. Previous dates of imagery used should include: 

1) one of the images covering the average date of 
maximum snow accumulation in an average snow year; 

2) at least one image from an early season snow date 

in the current snow season; and 3) if possible, another 
image from a date well into the snow season represent- 
ing average snow distribution for that date. As will 
be discussed later in Step 5, these previous snow season 
images are used to maximize the correlation between the 
satellite-derived snow water content indices and ground 


measurements of snow water content for any given 
sample unit on the forecast date in question. 

Step 4 ; Snow areal extent is estimated by ISU for the Landsat 
snow season forecasting date of interest, again using 
the procedure described in the Snow Areal Extent Pro- 
cedural Manual. The length of the delay that can be 
tolerated between image acquisition by the satellite 
and imagery receipt by the user will depend on the use 
of the snow water content estimate. A very short 
turn-around time (hours) will be required for weekly 
or sub-weekly water yield forecasts, but longer 
times usually will be satisfactory for monthly or 
seasonal forecasts. 

Step 5 : Snow areal extent data is next tranformed to snow 

water content index values by Landsat ISU. Snow 
water content is estimated from the following first 
order, time-specific model, designed to reflect 
physical snow depletion behavior in a given melting 
environment : 

J 

\ = r (M. .) (G ) K. where: 

j=l J J 


estimated Landsat snow water content index 
for image sample unit i, correlated to 
corresponding actual ground snow water 
content data; 

midpoint of the snow cover class (see 
Snow Areal Extent Procedural Manual) 
for image sample unit i on date j as 
obtained by interpretation of Landsat 
imagery; midpoint is expressed on a scale 
of 0.00 to 1.00; 

weight (related to snow depletion zone 
behavior characteristics) assigned (0.00 
to 1.00) to a past M. . according to the 
date of the current estimate; 

the number of times out of J^that sample 
unit i has greater than zero percent snow 
cover; and 

total number of snow season dates used, 
number of current snow season dates used. 
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The basic assumption of this simple model is that the 
greater the sum of areal snow extent over all dates 
used and the more often snow is present, the higher 
will be the expected snow water content. This model 
is most appropriate in mountainous environments. Each 
user may want to specify a more sophisticated and more 
physically realistic model for his own watershed (s) . 

Definition of the weighting factor, Gj, will be water- 
shed specific. Generally, weights should approach 1.00 
as the dates of imagery approach the date of forecast. 
The justification of the weighting relationship is 
that successively more recent climatic events have 
progressively more importance in determining the actual 
snow water content on the forecast date in question. 

To insure reasonably high correlation between xi and 
corresponding ground snow water content values, at 
least three snow season dates are usually considered 
(J ^3). Of these, one or two dates of Landsat imagery 
would be required during the early snow accumulation 
season. Occasionally, J may be only two, such as when 
the first date consists on an April 1st snow water map, 
based on the past year's Landsat data, and the second 
represents the current early snow season date in ques- 
tion. In all cases the sample unit grids on all dates 
must be in common register with respect to a base 
date grid location. Starting weight suggestions 
are: 0,25 for the previous year image that is repre- 

sentative of an average maximum snow accumulation date; 
0.50 for the early snow season date; 0.75 for dates 
occurring a month or two before the forecast date; 
and 1.00 for the forecast date and for all dates 
occurring within one or two months prior to the fore- 
cast date. 

Step 6 : All ISO's in the watershed(s) of interest should be 

classified into snow water content strata. Strata 
are defined by dividing the range of ISU snow water 
content index values for a mid- snow season forecast 
date (e.g., April 1) into two or three natural group- 
ings (seen by plotting snow water content index versus 
the number of ISU's). If no natural groupings are 
present, then the range of ISU snow water content 
index values should be divided into two or three 
ranges having approximately equal numbers of TSU's in 
each stratum. ISU's are assigned to the stratum 
having the snow water content index range that includes 
the snow water content index value for the particular 
ISU at the mid- snow season forecast date used. This 
ISU stratification is achieved by first applying 
Steps 1-5 to a previous snow season that was as close 
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to average as possible and for which good Lands at 
pass| coverage is ' available. 

After each ISU has been coded according to its stratum, 
the number of ground sample units fGSU's) consisting 
of snow courses required to satisfy the user precision 
criterion is calculated by stratum. The procedure used 
for this calculation is presented in Appendix I. 

j 

Initial implementation of this procedure requires infor- 
mation regarding 1) tabular summaries ^ of individual 
ISU snow water content index values; 2) the accumulated 
Landsat snow water content index in each stratum; 3) the 
ground-based snow water content variability by stratvim; 

4) the correlation between ISU and GSU data; 5) the 
total number of ISU's per stratum; and 6) average ISU 
and GSU costs by stratum. 

The ISU and accumulated ISU snow water content index 
information is obtained from the set of x^»s generated 
in the stratification process on the earlier snow 
season L^dsat data set. 

The total number of ISU's per stratum, %, will be 
known from the statification. An estimated correlation 
coefficient betwesen ground snow writer content estimates 
and Landsat-based snow water content index values is next 
calculated by obtaining ground snow; course data matched to 
a sraal’i. sample of ISU's. The snow course data must have been 
obtained within two or three weeks of the mid-season snow 
forecast date used to derive the corresponding Landsat ISU 
values. Ideally, separate correlation coefficients should be 
computed for each stratum if a sufficient number (e,g.,_10) 
of snow courses exist for each stratum. Otherwise, a common 
correlation coefficient (ignoring strata) should be computed 
by the usual formula (see Appendix I) and 
^plied to all strata. 

An estimate of snow water content variability cart^B^ 
obtained by applying the standard formulasfor simple 
random sampling variance to the snow course/ sensor 
data used to calculate the correlation coefficient. A 
separate variance should be computed for each stratum. 

If only a very small number (e.g., ^ 5) snow courses 
fall in a given stratum, then the variance among all 
Xi values in the particular stratum can be substituted 
instead. In this case, the values must be scaled 
to ground snow water content values by use of a ground- 
to-Landsat regression equation based on matched GSU 
and ISU data. 
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Image preparation and interpretation costs must be 
documented also for use in calculating sample size. 
Labor, materials, and overhead costs should be docu- 
mented on an ISU basis as in the example given in 
Table: I. Cost per GSU (snow course) measurement 
can be developed from the individual user organi- 
zation's pre-existing data. All cost data for sample 
allocation should represent operational costs. 

For tlhe first year of actual implementation, the 
calculated number of GSU's per snow water content 
index stratum should be allocated for given ISU's in 
those strata. This allocation is based on equal 
probability selection from the tabular summaries of 
ISU's generated in Step 6. Consequently, ISU's 
are numbered from 1 to Nj^ in a given stratum and 
a random number table used to select the calculated 
number of ISU's to which GSU's will be matched. 

The snow course or sensor should be located in a snow 
accumulation and depletion environment representative 
of the average of such environments covered by the 
given ISU. If it is not efficient to reallocate an 
already established network of snow courses and sensors 
then the snow water content stratum should be identi- 
fied for each such ground unit by determining the 
stratum associated with the ISU covering that ground 
location. After three to five years of using this 
remote sensing-aided technique for snow water content 
estimation, the user may find it desirable from a 
precision (i.e., variance control) standpoint to 
reallocate some ground courses to achieve the optimum 
stratum-by->*stratum GSU sample sizes calculated previous 
ly. However, before this is done, updated GSU sample 
sizes should be calculated using stratum- specific 
ground or image snow water content variance data 
averaged over successive seasons. 

Estimates of watershed or sub-watershed snow water 
content should be calculated, using the equation 
described in Appendix II. The values thus obtained 
should be entered into statistical or physical models 
to predict water yield. For example, the user could 
employ the remote sensing-aided snow water content 
estimate as an input variable in a regression equation 
for predicting water runoff. 

Finally, the utility of the above -described procedure 
for using remote sensing as an aid to estimating snow 


water content should be evaluated in regard to satisfy- 
ing the water yield forecasting organization's legal 
or contractual requirements. Utility can be judged 

by new information gained (e.g. , watershed physical 
relationships), forecasting accuracy or precision im- 
provement, cost savings, or forecast timeliness. In- 
formation gained on the preci,sion/cost effectiveness 
of the procedure and the sensitivity of the water 
yield forecast models to the resulting snow water con- 
tent estimates can be used to: 

1) refine the model (Step 5) used to calcu- 
late the snow water content index; 

2) better define snow water content index strata; 

3) generate better GSU sample sizes and alloca- 
tion strategies; and 

4) refine the actual snow areal extent inter- 
pretation, imagery enhancement, and analysis 
procedures used in Steps 2, 3, and 4. 
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Table I, Example calculation of image sample laiit costs for 
a Landsat -aided manual snow water content inventory! 


Total Cost Cost per ISU^ 

I. Pre- Inventory 

A. Image Acquisition 

3 LANDSAT dates with 
3 bands per date 
@ $3 per band; the costs 
of 2 of these dates amortized 

over 5 dates $12.60 $.006 

Resource Photography 

(Medium Scale .Aer^i'a! 

Photography for Image 
Analyst Environmental 
Type Training) 

B. Image Sample Unit 

Grided Landsat Color Composite 

Print Generation 

Film processing, and 

Printing > 

3 dates @ $11 per date 
The costs of 2 of these 

dates amortized over 5 dates $15.i»0 $.007 


Labor i 

0.5 hours per date |$13. 50/hr 
including overhead, 3 dates, 
the costs of 2 of which are 

amortized over 5 dates $ 9.^5 $.00-4 


1, Cost data based on 1975 University of California figures. 

2, Cost per image sampl‘d unit assuming 2218 (780, 000^^ test watershed) 
image sample units in the watershed (s) of interest. 

3, Two $500 flights amortized over 5 years, 7 dates per year, and two 
watersheds . 


$14.29 


$.006 


Table I, Continued 
II. Inventbry 


Total Cost 


Cost per ISU 


A. Interpreter Training 

1 hr per date, @ $13, 50/hr 
3 dates, the costs of 2 of 
which are amortized over 
5 dates 


$ 18.90 


$.009 


B. Image Interpretation 

Ave. 6 hrs per ‘date 
@ $13. 50/hr (2218 Image 
Sample Units) 

3 dates, the costs of 2 of 
which are amortized over 

5 dates 

C* Data Keypunching 

6 firs per date @ $13. 50/hr; 
3 dates,, the costs of 2 of 
which are amortized over 

5 dates 

D. Computer Analysis of 
Image Analyst Results 

0.075/hr @ $I»0/hr 

E. Selection of Random 
Numbers to Define Ground 
Sample Units 


0.5/hr @ $13. 50/hr amortized 
over 5 dates @ $13. 50/hr 


$ 113.^0 


$.051 


$ 113.^0 


$ 3.00 


TOTAL 


$1.35 


$301.79 


$.051 


$.001 


$.001 


$.136 
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Appendix I 


DETERMINATION OF GROl'ND SAMPLE SIZE 


The groimd sample size in each stratiim will be defined to be 
directly proportional to the fraction of total snow water content 
in that strattan relative to the total snow water volume over the 
entire watershed. This sample allocation method, known as X-pro- 
portional (see Raj 1968:65-68), uses the Landsat snow water content 
index data as the measure of size. Thus the number of ground 
snow courses required per stratum is given by 


h X ’ 


(equation 1) 


where 


h = the stratian index, 

n = the ground sample size over all strata, 

Xj^ = the sum of Landsat snow water content index 

values for all image sample units in stratum h 
and, 

X = the sum of Landsat snow water content index 
values for all image sample units over all 
strata. 

The total ground sample size, n, can be calculated by substituting 
the expression for nh given in equation 1 into the standard formula for 
variance of a stratified sample. Thus 


A A y ^ 0 ^ O 

X^'-yh -fi Vyh ' 


(equation 2) 


where 


A A 

V(Y) = the estimated variance of the final estimate 
(f) of total snow water content on the water- 
shed, 

Nj^ = the total number of Landsat image sample units 
in stratum h. 


< a 


n,X,Xj^ = as defined previously, and 
'>2 

a = the estimated variance for the 
^ estimate of snow water content in 
stratum h. Mathematically, 

2 


^ - ^'h {equation 3) 


where 


S , = the estimated variance among ground snow course 

^ measurements of snow water content expressed in 
inches of water, i.e.. 




where 


'ih 

^h 


average snow water content for snow course i in 
Stratum h, and 

average estimate of snow water content among 
snow courses in stratum h, ' 


p. = the estimated correlation coefficient between Landsat 
image sample unit snow water content index values (x-,) 
and spatially matched ground measurements of snow waier 


content in stratum h, i.e. 


h “ - V • - yh> 


where 


Xj^ = the estimated average Landsat snow water content 
index value per image sample unit in stratum h. 
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After rearranging terms in equation 2, the total grounci sample size 
can be expressed as 



(X T. 




Vyh ’ 


V(Y). 


(equation 4) 


Equation 4 is used to determine n and subsequently nj^ when the 
final estimate of variance (V(Y)j to be achieved is fixed at some 
desired level. Therefore, to solve for n, the user must first 
choose a value of ^(Y) in equation 4. The size of V will depend on 
the amount of sampling error in the final estimate of snow water 
content that can be tolerated by the user. For comparative pur- 
poses, the size of this error can be expressed as the ratio of 

J V(Y) to Y, that is, the standard error of the estimate of water- 
shed snow water content expressed as a fraction of the estimate. 

The smaller this ratio, the greater the precision of the final snow 
water content estimate. 


It should be noted here that the estimation process itself is 
assianed to be unbiased, given accurate ground and image sample unit 
measurements. That is, the expected value of the watershed snow water 
content estimate (obtained over all possible ground samples) is assumed 
to equal the actual value. Thus the sample size calculation only 
considers the error entering through the act of sampling a small 
fraction of the ground area. 


In addition to fixing V, the user must provide estimates for 
the variance in each stratum, These may be available from 

earlier studies. If not, they must be estimated by using equation 3. 
In this case, ^ is calculated from ground snow measurements taken 
in each stratum, pj^ using image data matched spatially with the ground 
data*, and guesses are made for n and n^. After solving for n (using 
equation 4) and then nj^ (using equation 1), the revised estimates of 
n and nh are substituted back into equation 3 and another solution is 
obtained for n and n^. This iterative process should be continued 
until the resulting estimates of n and n^ stabilize. 


If the user chooses not to specify a desired final variance V, 
but rather would like to minimize V subject to a limiting budget 
constraint, then the following formula for n should be used: 

n = ^ ~ (equation 5) 

X 


*If only a very small number of matched ground and Laudsat data pairs 
are available in each stratum, then these pairs will have to be 
pooled to calculate a common p among strata. 
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In equation 5, 


C - the total budget available for the snow water content 
survey, 

Cq = the total fixed (non -variable) costs in the inventory 

including the cost of generating snow water content index 
estimates for each image sample unit on the Landsat imagery 
within the watershed(s) of interest, 

C}j = the average cost associated with travel to and from, and 
measurements taken at a ground snow course in stratum h, 
and 

Xjj and X = the quantities defined in equation 1. 

Once X}j abd X are obtained from the Landsat imagery, equation 5 is 
used to solve for n and equation 1 is used to solve for the ground 
sample size, nj^, in each stratum. 
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Appendix I I 


FORMULAS FOR ESTIMATING WATERSHED 
SNOW WATER CONTENT 


A linear regression equation is used to combine ground and 
Landsat image sample unit information to produce a final watershed 
estimate of total snow water content. This estimation model is 
expressed as follows for a given stratum: 


’^h ’ \ \ = \ '*h * '’h (equation i) 


where : 


h = stratum index, 

A 

= total estimated snow water content Cexpressed as a 
volume of water) for stratum h, 

A 

Yj^ = estimated average snow water content (inches* of 
water) per image sample unit for stratum h, 

= area per image sample unit, 

Nj^ = total number of image sample units in stratum h, 

= estimated average snow water content (inches of 
water) per groinid sample linit for stratum h based on 
a sample of ground sample units, 

Xj^ = estimated average snow water content index per image 
sample unit for stratum h based on a sample of image 
sample units corresponding spatially to the ground 
units sampled, 

A 

Xj^ = estimated average snow water content index per image 
sample unit for stratum h, based on the estimates from 
all image sample units in strattim h, and 

A 

bj^ = estimated regression, coefficient between and Xj^. 


i.e., bj^ - " V ^^ih " ^h^^ 


% ^ : 
"f Cxih -x) 


where: 


Xih = the Landsat image snow water content index for image 
sample unit i in stratum h, 

or equivalent metric measure of depth. 
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^ih = the ground snow course snow water content (expressed 
in inches of Water) for the snow course located in 
image sample unit i of stratum h, and 

j\ j... A 

and are defined as above. 

The stratum-specific estimates may be combined to produce an 
estimate of total snow water content (expressed as a volume of water) 
for the ehtire watershed of interest by the following equation: 


Y* * 1 = S Y, 
total h 


(equation 2) 


where 


L = total number of strata. 
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Preface 

In this study, an attempt was made to distinguish areas of different 
surface environments from Landsat data, specifically in terms of vegetation, 
soil and terrain, and to evaluate the utility of Landsat data for such land 
surface classification. The factors of vegetation, soil and terrain (here 
meaning the combined effects of slope, aspect, and elevation) are primary 
in governing the amount of reflected electromagnitic radiation from a given 
area. Solar radiation incident on an area is differentially reflected due 
to the nature of the terrain, soil and vegetation; consequently, there tends 
to be a unique reflecteds'.signal that is characteristic of a specific com- 
bination of these factors. This study investigates the feasibility of 
differentiating land surface systems in a specific wildland area of approximate 
ly 41,250 acres (64. 5 square miles) within the Lassen County portion of the 
Lassen National Forest near Susanville, California. The information contained 
in this study is based on the analysis of reflected radiation sensed by the 
Landsat multispectral (MSS) imaging system. 

Classification 

When one deals with a task such as land surface investigation, there 
is added difficulty in classification due to the obscure definition of the 
individual components. iVhat might be identified as a different land form 
or surface type may be only a transitional form between larger, more cons- 
picuous or distinct forms. The problem of differentiating individuals can 
be somewhat reduced by looking at more than a single type of differentiating 
characteristic. According to R.L. Wright (1972), the principal requirements 
are that differentiating charactistics should be: 

1) Intrinsic properties of the things to be grouped; 

2) Observable, measurable, readily accessible, and ideally of a 
comparatively permanent nature; 

3) "Motivating" properties with the maximum number of accessory char- 
acteristics; 

4) Important for, and capable of application to, the objectives of 
the classification; and 

5) Suitable for the construction of a hierarchy. 


In addition, the classification scheme itself must be tailored 
closely to the purpose and objectives for which it was designed. The 
above basic considerations were instrumental in the design of this 
investigation. 

Approach 

In the following study a two-pronged approach was utilized in 
the analysis of land surface types. The analysis sought to discern 
areas of physical continuity on the basis of land surface appearance 
coupled with the aspects of vegetation and land use. 

The first part of the investigation dealt with the experimentation 
and use of various Landsat-1 MSS band combinations and permutations 
in the form of band ratio color composites as an aid to the manual 
interpretation of different land systems. The second part utilized 
a completely objective machine approach to land system differentiation 
and classification. A minimum of human interaction was used in the 
machine processing and classification of the study area because it 
was desirable to determine whether a purely automated approach to 
land classification produces results comparable to those obtained 
by the manual approach, or whether the automated approach is even 
applicable to such a problem. 

Discussion and Analysis - Part I 

A technique that combines digital computer processing and color com 
positing has been utilized to enhance the subtle spectral reflectivity 
differences often characteristic of land surfaces sensed at different 
wavelengths (Aviation Week and Space Technology, October 1977) . This 
procedure is called band ratioing. 

In the process of band ratioing, features which have little 
visible difference compared to their surroundings can be enhanced 
by emphasizing the relative responses between different spectral 
bands. Band ratioing eliminates common brightness trends between 
spectral bands and emphasizes spectral reflectivity differences 
sensed at different wavelengths. Normally different substrates have 
different reflectivities in the visible and near-infrared spectral 
regions. Differences are generally quite subtle and are difficult 
to discriminate on standard MSS images. In addition, the area being 
analyzed in this study possesses significant local relief. Band 
ratioing generally removes first order brightness caused by such 
topographic, slope and aspect variances (Rowan, et al., 1975). 

The actual process of band ratioing requires unique analysis 
equipment that can be found and utilized in only a very few locations 


across the United States. Due to the extremely flexible nature of the 
equipment of the Remote Sensing Research Program (RSRP) on the Berkeley 
Campus, such a technique was possible within the University of California 
system. 


In the ratioing process, the pixel reflectance sensed in one band 
is divided by the reflectance measured for that pixel in another band. 
This is done for each pixel. The resultant ratios, which may be 
displayed planimetrically in a color code, will show the variations 
between the slopes of the spectral reflectivity curves of the two 
wavelength bands. In some instances, differences in albedo are 
suppressed, and very dissimilar materials easily separable on a 
standard photographic image may be inseparable on a ratioed image 
because their spectral reflectivity curves in both wavelength bands 
are similar. Conversely, a distinct advantage of this method is that 
one type of material will appear the same or similar in a ratioed 
image regardless of the local topographic slope angle, since first 
order brightness trends are usually removed through band ratioing. 

Color compositing techniques offer an effective and visually 
pleasing means of combining ratioed images with standard MSS band 
images. Up to three bands may be combined utilizing blue, green and 
red filters. This is a particularly efficient method since two 
spectrally different areas may be nearly indistinguishable in a black- 
and-white ratioed image, while proper color combinations of up to 
three images permits discrimination on the basis of color differences. 
Manual discriminative ability is greatly increased because information 
from more than one image can be combined and because the human eye 
is capable of differentiating many more hue values than gray values. 


A number of ratios and ratio color composites were tried in 
dealing with the Lassen study are^; they included the following:* 


I) 

MSS 

7/MSS 

5 

Red 


MSS 

7 


— Green 


MSS 

6 


• — Blue 

III) 

MSS 

7/MSS 

4 

--- Red 


MSS 

6/MSS 

5 

— Green 


MSS 

5/MSS 

4 

Blue 

V) 

MSS 

7/MSS 

5 

-— Red 


MSS 

7/MSS 

4 

— Green 


MSS 

6/MSS 

5 

Blue 


II) MSS 5/MSS 4 — Red 
MSS 6/MSS 5 —Green 
MSS 7/MSS 4—- Blue 
IV) MSS 6/MSS 5 — Red 
MSS 7/MSS 5 — Green 
MSS 7 Blue 


(Figures X-1 through X-5) 


See pages 4 through 6 fot color ratio composite illustrations. 





Figure 


1. Band-ratio color composite #1 


MSS 7/MSS 5 
MSS 7 
MSS 6 


Red 

Green 

Blue 



Figure 2. Band-ratio color composite #2. 

MSS 5/MSS 4 --- Red 
MSS 6/MSS 5 --- Green 
MSS 7/MSS 4 -— Blue 
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3. Band-ratio color composite #3. 

MSS 7/MSS 4 — Red 

MSS 6/MSS 5 — Green 

MSS 5/MSS 4 — Blue 
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Figure 4. Band-ratio color composite #4. 

MSS 6/MSS 5 -— Red 
MSS 7/MSS 5 --- Green 
MSS 7 — Blue 
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Band-ratio color composite #5 

MSS 7/MSS 5 --- Red 
MSS 7/MSS 4 — Green 
MSS 6/MSS 5 — Blue 






Each band ratio color composite was analyzed using available 
geologic maps, high flight color infrared photography, natural color 
aerial photography, and on-the-site conventional ground photography. 

A discussion of each ratio composite and its individual components 
follows: 

Band ratio color composite #1 - MSS 7/MSS 5 Red 

MSS 7 — Green 

MSS 6 — - Blue (see page 4) 

I ' ■ , ; ^ i 

'In general, this particular combination of bands seemed to dis- 
criminate the largest number of surface types with the most meaning 
and accuracy. The MSS 7/MSS 5 component ratio of this composite is 
of particular interest since it was chosen as being representative 
of surface flora. 

The reasoning behind the use of th^ MSS 7/MSS 5 ratio is that, 
as seen from an aerial or space view, vegetation is almost entirely 
represented by its foliage (i.e., leaves). Other parts of the 
plant (stem, flowers, fruits, etc.) are usually not visible from 
above to the degree' that the leaf canopy is. Hence, in remote sensing 
it is the leaf that determines how an individual plant, plant asso- 
ciation, or community appears on a photo or image. 

Within the range of wavelengths (6.5 - l.ly) to which the Landsat 
MSS is sensitive, foliage typically interacts with solar radiation in 
the following manner: "Reflection or emission from the leaf’s cuticle 
and epidermis is relatively minor, and ordinarily need not be con- 
sidered; most of the light from the blue and red ends uf the visible 
spectrum is absorbed by chlorophyll while green light is largely 
reflected by it, thus accounting for the green appearance of the leaves 
Energy from the near infrared part of the spectrum is little affected 
by the chlorophyll, but is greatly affected by the gross structure of 
the spongy mesophyll tissues deep within the leaf. Just as snow and 
foam appear white because, with their spongy structure, they are 
good reflectors of energy, so the spongy mesophyll is also a good 
reflector. In fact this tissue appears very white when the over- 
lying portions of the leaf are dissected so as to expose it, but as 
long as the overlying tissues remain in place, the radiant energy 
which penetrates the mesophyll (and thereupon is reflected from it) 
is of far greater intensity in the near infrared than in the visible 
part of the spectrum." (Colwell, 1968) 

From the foregoing; statement it is obvious why MSS band 7 
(near infrared) was chosen as part of the ratio to represent the 
vegetation aspect of the environment. Differentiation of forest 
types is most easily accomplished using MSS band 7 where reflectance 
separation is great enough so that discrimination is possible. Also, 
the low values in MSS band 5, coupled with its fluctuation as related 


to plant vigor, indicate that a combination of MSS bands 7 and 5, 
that is, one ratioed with the other, v\rill give the image interpreter 
an idea of both vegetation type and vigor. Other research (Sabins, 
If'yS) has indicated that ratioihg two bands not only enhances inter- 
pretability in terms of the information represented in the two bands, 
but also aids in eliminating inconsistencies in the spectral data 
caused by differing slope, aspect and sun-angle relationships. 

The other two components of this first ratio composite are MSS 
bands 7 and 6. MSS band 7 is included simply as an additional factor 
in aiding vegetation discrimination. MSS band 6 is used as a spectral 
indicator of differences in the topographic elements of the landscape 
(i.e., slope and aspect) including soil. By registering and overlay- 
ing these three components, assigning each a primary color (red, green 
or blue), one can create a color composite that represents a pixel - 
by-pixel (grid method) classification of the landscape elements. The 
final product may be both manually (qualitatively) and quantitatively 
interpreted to determine the accuracy and meaning of such an infor- 
mation display method. > 

In the present test, through a close visual analysis of this 
first ratio composite, eight distinct categories of land surface 
types (or photomorphic units) could be discerned (see Table 1) . 

Each of the eight land surface types defined by its color charac- 
teristics on the ratio composite was quantitatively analyzed to 
determine its spectral characteristics. As seen from the statis- 
tical summaries of each of the areas (see Tables 2 through 9), 
there is a histogram of pixel reflectance values for each band in 
each of the eight categories. From each of the histograms a mean 
(representative) reflectance value can be derived that will represent 
the reflectance value for that land surface type in that particular 
band or band ratio. 

In summary, each of these areas is unique and represents distinct 
elements of the landscape. Virtually the entire study area may be 
described through these eight basic land systems. Deviation from 
any one of the eight basic units seems to occur either in transition 
zones between units or in areas that are of such great topographic 
irregularity that reflected radiation does not represent the true 
surface character (e.g.. Fault Zone, Antelope Mountain). Although 
such a segregation of land units seems to work reasonably well, it 
must be reiterated that the specificity of the classification is 
dependent upon the classification’s designated puipose. For some 
purposes an eight-element land classification scheme such as this 
would be perfectly adequate. For instance, a type of large-scale, 
regional land surface inventory may be accomplished quite well by 
using the eight basic land units to build on a higher order classi- 
fication scheme. Obviously detailed vegetation/habitat discrimination 
is impractical using such a coarse scheme, but it is conceivalde that 
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Table ,1 

LAND SURFACE CATEGORIES 
Ratio Composite 1 


ORIGINAL PAGE IB 
OF POOR OITAtU^' 


AREA 

TONAL DESCRIP. 

GEOLOGY 

SLOPE/ASPECT 

VEGETATION 

1 

Blactc 

Qru^’-Black 19.97fi/N-NW 

vesicular 
scoriaceous 
flows of divine 
basalt, charac- 
terized by lack 
of erosion or 
weathering, flow 
surfaces and 
structures locally 
well preserved 

Generally non- 
vegetated 

2 

Yellow, red 
brown mixture 

Qru*^ 

10.755i/N-NE 

Yellow Pine 
substrate 
not obscured 

3 

Red brown- 
almost exclu- 
sively 

Qru^ 

9.98%/E-NE 

Yellow Pine 
shrub, underbrush 
substrate obscured 

4 

Generally 
white with 
red brown, 
yellow and 
blue anamolies 

Qru^ 

15.36%/E 

Much less dense 
Yellow Pine forest 
with substrate 
obscured 

5 

Green brown, 
yellow, some 
blue 

Qru^ 

9.22%/E 

Much less dense 
forest with sub- 
strate evident 
through canopy 

6 

Intense blue- 
homogeneous 

Qpv*^-Pyro- 

clastic 

basaltic 

1 . 54%/S-SW 

Grasses/low sage, 
some pine, spruce 


ivudi ijr 

agglutinated, 
vent tuffs, 
poorly consoli- 
dated rhyolite 
pumiceous tuff 


9 
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Table 1 (continued) 

• "V 'l ; 

. V* * 

AREA 

TONAL DESCRIP. 

GEOLOGY SLOPE?' ASPECT 

VEGETATION 

t 

; .. ' 7 

i 

Light blue- 
homogeneous 

Qal-lpecent 0.77K/S 

alluvium 

Low sage/grasses 

! ' ' 


alluvial silt, 

M 

i 

i ' 

f . 

1 

i 

j 


sand, and locally 
coarse gravel, 
deltaic slope- 
wash stream 
channel, flood- 
plain deposits, 
fans partially 
lacustrine bog 
and swamp deposits. 



glacial outwash 




8 

Light blue with Qpv 

9.22?i/NW 

Shelterwood cut 


yellow and red 


Yellow Pine - 


brown interspersed 


recovery 

9 

Black (water) None 

None 

None 
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Table 2 

RATIO COMPOSITE 1 - FIELD 1 
Summary Statistics 


ORIGINAL PAGE IS 
OF POOR QUALITY 
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BAD DATA 
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Table 2 (continued) 


♦From this type of normal distribution, one would postulate that there 
is either one type of vegetation or n£ vegetation in this land unit that 
dominates the character of the reflected spectral signal. An examination 
of the band 7 reflectance histogram reinforces this idea. A further 
examination of high-flight photography and on-the-ground reconnaissance 
indicates that this area is largely void of any vegetation due to its 
physical make up and relatively youthful geologic state (basaltic flow). 



Table 3 (continued) 




f ■■ , 'Af 

♦This broken Gaussian distribution here seems to indicate that the surface 
is either covered by a number of co-dominant plant species/associations or 
that the substrate is greatly influencing the spectral signal. Examination 
of ground truth data substantiates the latter, since this area is dominated 
largely by one vegetation type that is not dense, and so the underlying 
substrate plays an important role in the character of the spectral signal 
reflected from this land unit. 
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Table 4 

RATIO COMPOSITE 1 - FIELD 3 

Summary Statistics OKIGINAL PAGE IS 

^ OF POOR QUALITY 



BAND 

MEAN 

STDEV 

CV 

SKEW 

•'KURTOSIS 

BAD DATA 

Blue! 

3(7/5) 

31.45 
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Red 
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Table 4 (Continued) 
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Table 4, Concluded 


♦The bimodal distribution presented in the band 5 histogram indicates the 
influence of two 'dominant factors on the spectral signal. An examination of 
the site through an on-the-ground visitation showed the two dominant features 
to be a yellow pine and a big sagebrush community. 
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RATIO COMPOSITE 1 - FIELD 4 
Summary Statistics 
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Table 6 

RATIO COMPOSITE 1 - FIELD 5 
Summary Statistics 
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Table 6, Concluded | 
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*The represented land unit system illustrates a high dominance by one *»sf | 

factor of the environment. Field investigation shows this area possesses | 

a sparse Yfellow Pine forest cover with little or no underbrush, and a ■»T | 

high percentage of the substrate exposed; hence the nearly Gaussian dis- , | 

tribution. I 
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Table 8 

RATIO COMPOSITE 1 - FIELD 7 
Summary Statistics 
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Table 8, Concluded 


♦This type of bimodal distribution would seem to indicate a dominance 
by two vegetation complexes or individuals within this land unit. Field 
reconnaissance has shown this to be the situation. The land unit is 
a mosaic of areas of exclusive grass dominance and areas of sagebrush/ 
grass association dominance. 


Table cj 

RATIO COMPOSITE 1 - FIELD 8 
Summary Statistics 



BAND 

MEAN 

STDEV 

CV 

SKEW 

KURTOSIS 

BAD DATA 

Blue 

3(7/5) 

36.89 

1.54 

4.18 

-.24 

-.23 

0 

Green 

4(7) 

21.30 

.93 

4.37 


-.43 

0 

Red 

5(6) 

.92 

.09 

9.48 

.43 

2.32 
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I 



Table 9> Concluded 


♦The nature of this thinned (shelterwood cut) area is one of a conglomeration 
of ground covering elements. Yellow Pine covers this land unit in a sparse 
distribution, while an understory of sagebrush with grasses covers part of 
the space in between, with yet a large percentage of the area being non-vege- 
tated. This probably accounts for the normal yet slightly dispersed distribution 
of reflection values. 
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a method of band ratio compositing could be developed for such a 
purpose. Again, this elementary land system breakdown was done only 
for I the area designated in this study, and therefore, this is the ^ 
only area for which it can be used with confidence. By gambling to 
a limited degree, this eight -unit designation scheme may be extra- 
polated to an expansive area adjacent to the test site that is of the 
same general climatic/topographic make-up. 

Band ratio color composite #2 - MSS 5/MSS 4 ^ — Red 

MSS 6/MSS 5 — Green 

MSS 7/MSS 4 Blue (see page 4) 

Tliis particular composite represents a conglomeration of ratios 
that involve combinations that have been used in earlier geologic 
investigations (MSS 5/MSS 4, MSS 6/MSS 5) (Short and Marrs, 1975), 
as well as a ratio (MSS 7/MSS 4) that was developed in this study. 

The! band 7 to band 4 ratio was designed to illustrate a maximum 
amount of contrast between vegetated and unvegetated surfaces, 
while the other two components are essentially representative of 
varying substrate composition. Within the band 7 to band 4 ratio, 
band 7, as noted previously, is a good representative of living vege- 
tation because of the high infrared reflectance of live vegetation. 
Band 4 represents a secondary peak in vegetation radiation reflection. 
Hence, in the presence of vegetation, both the band 7 and band 4 
peaks will be present and result in a relatively high 7/4 ratio. 
Conversely, a lack of vegetation will result in a much lower band 
reflectance, and a lower band 4 reflectance (but not to the same 
degree as band 7) resulting in a lower 7/4 ratio. 

As shown on page 4, the riparian vegetation that dominates the 
stream channel directly south of Crater Mountain is displayed with ' 
a high 7/4 ratio (coded blue in this illustration). Relatively 
sparse vegetation on the flanks of Crater Mountain itself shows 
a light blue to green cast signifying the presence of the vegetation 
canopy along with a substantial part of the substrate coming through 
to influence the signature. Probably the most striking aspect of 
this composite is the way in which the band 5 to band 4 ratio picks 
up exposed geology. Along the upper flanks of Crater Mountain the 
exposed basaltic flows are given a reddish cast, while the lower 
elevation alluvial flats, made up of alluvium and exposed volcanic 
rock, show up as an intense red. The incompletely regenerated timber 
cut (on Antelope Mountain) registers a reddish cast showing the sub- 
stantial influence of the substrate on the spectral signature. The 
overall whitish cast given to a majority of the vegetated areas 
(especially on the flanks of the surrounding mountains) indicates 
an approximate equal contribution of color by each of the band ratios 
in the composite. This would seem to indicate that these areas are 
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vegetated, but the substrate is playing at least an equal role in 
affecting the reflected spectral signature. 

■ . 1 . 

Band ratio color composite #3 - MSS 7/MSS 4 --- Red 

MSS 6/MSS 5 --- Green 

MSS 5/MSS 4 Blue (see page 5) 

This composite differs from the previous one only in that each 
of the band ratios was assigned a different color. As can be seen, 
the results are the same (i.e. , areas of corresponding color on ratio 
composite #2 also correspond in color on ratio #3, but are portrayed 
in colors different from those in ratio composite #2). The intention 
here was to show that although both composites present the same data, 
different modes of presentation may permit easier qualitative inter- 
pretation. For example, in ratio composite #2, the riparian vegetation 
is slightly easier to distinguish from the surrounding alluvial flat 
merely because of the adjacent tonal contrast. This suggests that 
there is some utility in experimentation with color schemes attached 
to each band ratio composite. To make as many distinctions as 
possible, all possible color combinations should be viewed at least 
once. 

Band ratio color composite #4 - MSS 6/MSS 5 — Red 

MSS 7/MSS 5 — Green 

MSS 7 Blue (see page 5) 

In this particular display! MSS band 7 and the 7/5 ratio again 
reinforced each other in representing the ground vegetation cover. 

All greens, yellows, whites, and light blues represent some type of 
vegetation cover. Darker blue correlates well with exposed substrate 
and volcanic rock and illustrates how MSS band 7 represents more than 
just vegetation, and should not be used solely as an indicator of 
vegetation. 

Band ratio color composite #5 - MSS 7/MSS 5 Red 

MSS 7/MSS 4 — - Green 

MSS 6/MSS 5 Blue (see page 6) 

As indicated by the band 7 to band 5 ratio, vegetation ; on the 
slopes of Crater Mountain is of a sparse nature, with its spectral 
signature heavily influenced by the substrate (as indicated by the 
6/5 ratio) . Vegetation is shown with a blue to blue-violet color 
while other areas are white showing the mixed influence of vegetation, 
substrate and topography. Exposed geology again seems to be highly 
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(^correlated to the 6/5 ratio, exhibiting a green color on this com- 
j; ^osite. 

In concluding this section, we should restate that the purpose 
in experimenting with different band r?tio combinations was to deter- 
mine if individual land units or systems could be distinguished, and 
to look more closely at the problem areas that needed to -be field 
checked. 

I Within each particular land unit, basaltic rocks were easily 
distinguished, exhibiting a dark tonal signature on all of the com- 
posites investigated. Alluvial areas were easily discriminated, as 
were all areas of exposed substrate. Vegetated areas were rather 
easily noted, but confusion with the substrate was encountered when 
the[ surface was not entirely covered by a vegetation canopy. Re- 
ferring to ratio composite #1, spectral histograms and plots pro- 
vided an added dimension in aiding the manual interpretation of the 
different ratio composites. Often, the 7/5 band ratio could be 
compared with band 7 to verify the presence, extent and diversity of 
vegetation. MSS band 6 was somewhat useful in illustrating substrate 
consistency. In the other sets of composites, ratios such as 6/5 and 
5/4 aided in surface geology interpretation, and 7/5 and 7/4 ratios 
were of considerable import in interpreting vegetation. 

Although the land units identified in this part of the study 
may have been on a somewhat small scale, they have clearly served 
to indicate the potential value and importance of more detailed cost 
intensive investigations along this same line. The amount of data 
displayed within a color ratio composite may never be fully recognized, 
and for this reason, additional experiments should be undertaken to 
exploit these information sources to the limit that existing technology 
permits. 

Another major channel of application involves the use of com- 
puterized classification routines that classify raw Landsat data 
with a minimum of human effort and interaction. The second part 
of this study will institute this method and compare the results 
with those obtained through the band ratio method. 

Discussion and Analysis - Part II 

The following analysis involves the use of two classification 
procedm'es in ordering into discrete spectral classes the spectral 
reflectance values of each of the four MSS bands for each pixel in 
the study area. The computer classification includes an unsupervised 
clustering algorithm (ISOCLAS) and a supervised classification 
procedure using a maximum likelihood decision rule (GALSCAN) . 

ISOCLAS (iterative ^elf-Organizing Clustering Algorithm) is a 
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large computer program developed at the Johnson Space Center (as 
IS0CL$3 and later adapted to the CDC 7600 system at the Lawrence 
Berkeley Laboratory by the personell of the Remote Sensing Research 
Program. Several controlling inputs allow the user to control the 
size and number of , Clusters developed by the program. (A cluster 
in this case represents a nxanber of pixels grouped together into a 
spectrally homogeneous group.) 

Aside from this minimal number of human inputs, the program 
acts on its own to classify the digital spectral information pro- 
vided to it. Because the program produces a classification -type 
cluster map on which each cluster requires post-processing identi- 
fication and interpretation, it is frequently referred to as an 
unsupervised classification scheme. 

In the case of Wxldland areas, such as the study area involved 
here, ISOCLAS is commonly run on a smaller area (sub-area) within 
the larger one. The sub-area, in the mind of the intei'preter, 
must possess a good cross section of spectral classes represented 
within the entire larger study area. The final outcome is that 
the small sub -are a has been classified into a number of clusters, 
each representing a set of linique reflectance characteristics. 

This information is then used to ''seed” or supervise the supervised 
classification algorithm, CALSCAN. 


CALSCAN, referred to as the maximum likelihood classification 
algorithm, requires some type of initial input data on which to base 
its final classification. This can be either in the form of a train- 
ing set (i.e., spatial sample (s) of digital data of a known ground 
feature selected by the user from which the spectral characteristics 
are calculated and used as supervisory data on which the classifi-. 
cation is based) or previously calculated spectral statistics, such 
as output from ISOCLAS. The input data is then used as a basis in 
the calculation of the maximum likelihood ratio within the algorithm. 
Briefly, the maximum likelihood ratio is a probability decision rule 
used to classify an area of interest (target) from multispectral 
data. In its simplest form, the likelihood ratio is Pt/Pj^, or the 
probability P of an unknown spectral measurement being classified 
as the target t, to the probability of an unkno\m spectral measure- 
ment being classified as background b. When Pt/Pb_^l, the decision 
is t; when Pt/Pb^l, the decision is b. Probability density functions 
are calculated from the spectral samples provided as input (training 
sets or statistics). As the number of training sets increases, the 
mathematical computation of the maximum likelihood ratio increases 
in complexity. As a result, computer analysis is necessary. The 
entire process, from the generation of training sets or statistics 
through the maximum likelihood classification, is known as automatic 
data processing (ADP) or automatic spectral pattern recognition of 
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multispectral remotely sensed data. 

Within this study, the order of events was much the same as the 
procedure just described. Initially, Crater Mountain was chosen as 
the test area (sub-area) that would be processed through the imsuper- 
vised classification algorithm (ISOCLAS) . The Crater Mountain area 
(200 pixels by 150 pixels) was chosen primarily because of its sub- 
Sjtantial local relief (as great as any feature within the entire 
study area), variety of slope and aspect , classes , and diversity of j 
vegetation and surface types which seemed to cover the entire spec- i 
trum of land classes possible for the greater study area. The out- 
put, in the fonn of a class map, depicts eleven different spectral 
groupings. Gross patterns and trends were recognized and differen- 
tiated, but no substantial effort was undertaken to interpret the 
exact significance of each spectral class at this point in the 
analysis. The statistical information that resulted from this 
classification was used as input into the supervised classifier in 
the form of means, standard deviation and covariance, and correla- 
tion matrices. The supervised classification routine (CALSCAN) 
was then run using these training statistics. Final output from 
the algorithm was in the form of a class map designating which class 
each pixel had been assigned to and a probability map siting 
the probability of misclassification for each pixel. The latter 
output product is currently unavailable, but is particularly useful 
for detailed inventory analyses in which it is important to know the 
probability associated with each class when estimating the areal 
extent dr physical volume of particular features. The output class 
map is shown on the following page (Figure 6 and a description and 
analysis of the results follows. 

I A multistage analysis procedure was utilized to determine the 
significance of each of the eleven spectral classes produced by the 
maximum likelihood classifier (CALSCAN) . "Multistage" in this case 
refers to the use of different scales of imagery and high flight 
photography to sample and verify areas of interest. In this part of 
the investigation, a standard Landsat color composite (scale 1:100,000) 
was utilized along with a 1:435,750 scale enlargement of a Landsat 
color composite, high flight RB-57 1:130,000 scale natural color 
and color infrared photography, and ground reconnaissance photos. As 
a supplement to the different types of imagery used, ancillary infor- 
mation in the form of geologic and topographic maps was also referred 
to in evaluating the eleven spectral classes. Table 10 shows the 
results obtained. 

As with the color ratio composite method of land surface iden- | 
tification, the validity of computer classification has been shown. 

The land surface types identified in Table 10 correspond closely to 
those land surfaces described in Table 1. Hence, for the type of 
small scale land surface phenomena investigated in this study, both 
methods appear to be comparable. Ideally, a combination of both 
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Figure 


6. CALSCAN classification summary. (Refer to detailed caption 
on following page.) 



gure 6. Caption, continued 
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Table 10 

LAND SURFACE CLASSES - AUTOMATIC CLASSIFICATION 


Class No . 
1 




5 



9 

10 

ri 


Description/Location 

Intermittent patches dis- 
tributed throughout the 
extent of the alluvial flats 

Generally confined to nor- 
thern slopes, and slightly 
higher altitudes (related 
to aspect and altitude) 

Located along stream 
channels within the 
alluvial flats 

Generally located in the 
transition zones between 
higher areas and the 
alluvial flats 

j 

Generally found at the 
interface of two areas 
dominated by a single class 
(i.e, 2, 6; 2; 4; 2, 7; etc. 


Found in higher areas, 
largely in conjunction 
with class 2 

Distributed throughout 
the alluvial flat areas, 
and exclusively confined 
to them 

Located almost exclusively 
on SE facing slopes 


Found on the upper reaches 
of Crater Mtn., in exposed 
rift areas 

Found in proximity with 
class 9 


Crater Lake (on Cr. Mtn.) 
Poison Lake (\l. of Cr. Mtn.) 


Identification 

Very rocky, largely i 
barren basaltic area 


Oensley forested. Pine, 
underbrush present 


Wet meadow and riparian 
vegetation 

Vegetational trans- 
ition zone (big sage/ 
bunchgrass) to forest 

Transitional zone 
exhibiting the effects 
of a location adjacent 
) to two highly dominant 
cover classes 
! 

Sparse/low density 
forest; underbrush 
largely lacking 

Alluvial areas covered 
with bunchgrass and low 
sage 


Confusion class, illustra 
ting the affects of a 
relatively low sun angle 
or a high relief area 

Unvegetated basalt flows 


Exposed basalt vegetated 
(underbrush only) 

Water 










methods (ISOCLAS/CALSCAN classification of the band ratio data) and I 

a qualitative as well as quantitative analysis of both products would | 

serve as the best method in a more detailed land surface type invest!- | 

gation. I 

Conclus ion | 

~ h" ' o i 

I 

1. Comparability and Use | 

f. 

The degree of validity and usefulness of the color ratio com- I 

posite method of land surface identification/classification and the | 

automatic data processing approach appear to be comparable. Each I 

can be used independently of the other if the facilities are not I 

available to generate both products, but as a method of checking 1 

and verification the greatest accuracy can be expected from using 

the two products in conjunction with each other. | 

2. ^plicability of the Techniques I 

The type of information generated here through the classification I 

of spectral data may be used as input into data base systems and as 1 

inputs into regional planning decisions ox' resource inventories. 1 

Models of the physical environment that require information on sur- 
face types (e.g., climatic process models, geologic/geomorphic ’ 

process; models, or water resource models used in predicting runoff) r 

could use this type of data. 

3. Further Investigation 

There is much opportunity for new and more in depth invest!- i 

gation concerning band ratio analysis and computer processing of 
spectral data. With adequate funding, a combination of the two 
methods involving experimentation with new ban.d ratios and later 
ISOCLAS/CALSCAN computer processing of the results could provide 
new insight into land surface differentiation. More detailed (larger 

scale) differentiation of land surface types could open new areas of I 

possible application of the technique. Included among these could I 

be soil group/association identification, geologic small scale 

mapping, detailed vegetative mapping, terrain analysis, water ; 

quality analysis, and the monitoring of disturbance activities. | 

(e.g,, mining and logging). | 
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